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ABSTRACT 


A procedure for the application of the Finite Element 
technique to analyze the coupled stress, mass and heat balance 
euations in a two-dimensional Orthoeonal metal cutting situation 
is presented. The material undergoing deformation is treated as a 
strain-rate sensitive, non-strain hardening visco-plastic solid. 
The visco-plastic behaviour is utilized to describe the plastic 
flow as equivalent to that of an incompressible, non-newtonian 
fluid. For the sake of simplicity, thermal softening is not taken 
into account. The resulting mass, momentum and heat transfer 
equations are solved iteratively with successive under-relaxation, 
using the Finite Element Method. The temperature and strain-rate 
distributions thus obtained are used in determining the shear plane 
and the dimensions of primary and secondary deformation zones. A 
parametric study has also been conducted to examine the effects of 
process variables such as the cutting velocity, chip thickness 


ratio and depth 

of 

cut on 

the 

temperatures at 

the chip-tool 

interface. This 

is 

because 

the 

high temperature 

distribution 


occurring in this region due to frictional heat dissipation is the 
prime reason for tool wear. 
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CHAPTER 1 


INTRODUCTION 


1 . 1 GENERAL BACKGROUND 

The Importance o£ machining in the modern engineering 
industry need not be over-emphasised when it is stated that nearly 
half of the engineering products are produced through machining at 
one stage or the other. In this age of wide-spread automation, ve 
are witnessing the trend that the demand for productivity and 
process optimization is Increasing day by day. As the search for 
new materials and processes is continuing, effort is still 
underway to improve conventional processes by reducing the 
material wastage to a bare minimum, increasing tool life and 
improving the product finish. In view of the crucial role played 
by machining among all the material-processing operations, a 
scientific understanding of various metal-cutting processes is 
vital for modern engineering technology and practice. 

A detailed analysis of the mechanics of metal-cutting which 
sheds light on the underlying plastic deformation processes and 
stresses la very helpful in the prediction of some of the input 
parameters such as the cutting forces. The evaluation of the 
cutting forces, in turn, enables one to determine the power input 
and to choose the tool material. Further, minimization of the 
power input for a certain material removal rate QMRR) often forms 
the objective of optimizing the metal -cutting operation. Also, an 



insight into the effects of various process variables like the 
tool geometry, tool material properties, cutting velocity and 
tool-tip temperature is useful in determining the optimal cutting 
conditions, for lowering machining costa and increasing the 
productivity. 

Over the past few decades, a large volume of experimental 
data has been collected on various aspects of metal-cutting, but a 
thorough understanding of this highly complex process is yet to 
emerge. Till date, a comprehensive theoretical treatment is also 
not available, though many an attempt has been made in the past in 
that direction. The present work, in its right earnest, is a small 
effort to fill some gaps existing in the theory of metal cutting. 
A coupled analysis of the deformation and thermal processes in the 
vicinity of the tool-tip has been attempted here, which is 
expected to provide an accurate estimate of the temperature field 
as well as the size of the plastic deformation zone. 

1.2 REVIEW OF PREVIOUS WORK 

Research in the area of metal cutting commenced right from 
the post-war period. Development both on the theoretical and 
experimental fronts took place gradually, with better theoretical 
models replacing their predecessors, as the computing and 
experimenting techniques advanced. 

The first attempts towards developing an understanding of 
the Mechanics of Metal Cutting were performed by analysing the 
chip-formation process. As early as 1945, a pioneering work by 
Mer chan tCl] and Pliapanen L2 ] , was presented on the so-called 



3 


classical "single shear plane model". This model assumes that the 
chip is formed due to plastic deformation occuring along a single 
plane, known as the shear plane. The most advantageous aspect of 
this model is that it enables the determination of the average 
yield shear stress and the slip velocity along the .shear plane, 
purely through simple geometric constructions. It also established 
beyond any doubt that metal cutting is basically a shear 
deformation process. The Merchant Piispanen model was 
subsequently improved by Palmer & Oxley [3] and Okushima & Hitomi 
[4]. These studies suggested that though the formation of the chip 
occurs due to the flow of the metal under shear, the deformation 
is not limited to a single plane; it occurs in a narrow region 
approximated by two parallel planes. This model later came to be 
known as the Thick Shear Zone Model. In view of the small 
thickness of the shear zone in comparison to its length, this 
model assumes that the state of stress within the deformation 
region is uniform simple shear. 


Around the 

same period 

of 

the 

above-mentioned 

theoretical 

models efforts. 

were going 

on 

to 

val idat e 

these 

models with 

experimentation . 

Kececioglu 

[5] 

was 

the first 

person 

to observe 


the process through photo-micrographs in 1958. The plastic zone 
where the chip formation takes place, was photographed at short 
intervals using a quick-stop device. By noting the deformation of 
the grain boundaries, the size and the shape of the plastic 
deformation zone were estimated. This study established that the 
plastic region could be approximately represented as a thin 
parallel-sided zone. Later in 1969, a more effective way of using 
the photo-micrograph technique was proposed by Stevenson and Oxley 



[ 6 ] . Asauming a thin parallel-sided zone, they observed printed 
grids on the material, before and after the deformation. From the 
streamlines of metal flow, the strains and strain rates were 
calculated. The important process variable of flow stress (called 
by different names such as effective yield stress, dynamic shear 
stress (DSS) etc.) could also be evaluated, as a function of 
strain by this method, for a wide range of strain rates. Later 
investigators [7,8] proposed that the flow stress is a strong 
function of strain rate as well. 

A totally different approach was proposed by J.T. Black [7], 
based on the theory of dislocation mechanics. It suggested, that 
the flow stress, is composed of two parts i.e. thermal and 
non-thermal . The observations of the study indicated that the 
plastic zone is divided into alternate shear and lamella bands. 
The shear bands (or shear fronts) give rise to the dominant 
non-thermal part while the lamella regions, lead to the thermal 
part of the flow stress. In this study, the flow stress was 
evaluated by relating the stress to the Stacking Fault Energy 
(SFE) , which itself was obtained from the dislocation distribution 
in the metal. This approach, mostly used by metallurgists, has not 
progressed far and and it is difficult to apply this approach as 
long as the concepts of dislocation mechanics do not match exactly 
with those of continuum mechanics. 

In all studies discussed above, the emphasis was only on the 
primary deformation zone occuring across the shear plane. It was 
Von Turkovich [8,9] who first observed the peculiar form of 
plastic flow, extremely localised in a small region in the chip. 
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ahead of the tool ed£e along the rake face. This was called the 
Secondary Shear Deformation Zone. The experimental studies 
conducted so far, have established that material behaves like a 
highly viscous fluid within the deformation zones and as a perfect 
solid (with infinite viscosity and zero strain rate) outside 
them. The viscosity, in such a fluid like statein fact, depends 
strongly on the material deformation rate. 

With the advent of high-speed computers coupled with the 
increased awareness of the immense potential of numerical 
techniques like FDM and FEM considerable interest was daoted to the 
application of these procedures. The credit for the effective 
usage of the FEU technique to visco-plastic analysis goes to 
Zi enki ewicz et al .£.12# 1 3} - In 1973 , these researchers applied FEM 
for the modelling of metal forming and extrusion problems. In 
metal cutting, the FEM technique was employed by Tay and Stevenson 
[14l^for predicting the temperature field for orthogonal-cutting. 
They used hyperbolic streamlines in the primary deformation zone 
for calculating the velocity fields and strain-rates and derived 
semi-empirical relations for computing the heat generation. The 
only drawback of this study was that it did not provide a coupled 
analysis of the plastic deformation and thermal processes during 
metal-cutt ing . 

In study similar to that of Tay & Stevenson[l4,3^, exhaustive 
experimentation backed by theoretical validation was conducted by 
Muraka & Hindu ja[ 17]. These authors obtained the average chip-tool 
interfacial temperature for a wide range of cutting conditions. 
Balaji [18] applied the FEM technique to evaluate the temperature 
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field within the work-material and the tool, for coated carbide 
tools. He also measured the average tool-tip temperature by the 
tool-chip thermocouple method and found reasonable agreement with 
the FEn predictions. 

The FEM solutions available at present for the metal-cutting 
situation provide only the temperature distribution in the chip 
and the tool for a given velocity field. In these analyses, the 
velocity fields themselves are obtained by empirical and 

geometrical considerations. The relations obtained from such 
considerations, however, suffer from gross simplifications which 
are necessitated by the complexities of the process as well as the 
domain. These simplifications, in turn, tend to alter the problem 
being solved considerably and also affect the accuracy of the 
solution. Also, the available temperature predictions through FEM 
have not been adequately compared with experimental results due to 
limitations in measuring the interfacial temperatures accurately. 
1.3 OBJECTIVES «. SCOPE OF PRESEMT STUDY 

In the present study, therefore, an attempt has been made to 
fill the gap in the analysis of metal-cutting by solving the 
coupled flow and energy equations. Stresses in the cutting zone 
have been modelled by considering visco-plastic material behaviour 
during deformation. Accordingly, the metal flow beyond the 
elastic limit is treated similar to that of a liquid of very high 
viscosity. The viscosity itself is described as a strong, 
non-linear function of the strain-rate and yield stress of the 
work-material. The material is considered to be purely strain-rate 
sensitive, with no effects of work-hardening or thermal softening. 
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A wore ceneral form of material deformation behaviour was not 
taken into account for the sake of simplicity. 

It is well established from previous experiments that the 
deformation is restricted to a very small region around the tool 
tip. Therefore, a solution domain extending upto a few millimetres 
from the tool-tip into the workpiece and the chip has been 
considered. Beyond this domain, the material is assumed to be 
perfectly rigid. The tool is also taken to be perfectly rigid and 
sharp with no built-up edge. The tool edge is assumed to be 
orthogonal to the cutting direction. Further, for generating the 
domain and the finite element mesh easily, the geometry of the 
problem has been simplified in the following manner. The chip is 
considered to be bounded by plane parallel surfaces beyond the 
contact point. Adjacent to the shear plane, the free-surface of 
the chip is assumed to be part of a cylindrical surface with 
constant radius of curvature. Though the chip-curling process 
makes the real boundaries of the problem very complex, the 
simplifications invoked in the present study are appropriate, as 
the focus of this investigation is to perform coupled stress and 
heat transfer analysis for a given domain. 

Another important assumption invoked in the present study is 
that all the power expended in plastic deformation is converted 
into heat. Though some residual str esses/strains are retained both 
in the workpiece and the chip, this assumption simplifies 
computation to a great extent without sacrificing the accuracy 
of predictions, when compared with the experimental results. 



The fioverninc equations of flow and enercy resultinc after 
the incorporation of the above-mentioned assumptions, have been 
solved usinfi the Finite Element Technique. The velocity, 
temperature and strain-rate distributions have been obtained 
and utilised to identify the shapes of isotherms and the plastic 


zone boundaries. 



CHAPTER 2 


THEORETICAL FORMULATION 


In many ways, metal cutting is unique among plastic 
deformation processes. The geometry of the deforming material is 
unbounded in metal cutting, which is not the case for most metal 
forming operations. Also extremely localised asymmetric 
deformation occurs at exceedingly high strains and strain rates. 
The metal cutting action itself is a consequence of high 
compression followed by shearing of the metal. Not to mention the 
least, the fact that all the deformation phenomena occur in a 
minute volume of a few millimetres size, adds the new dimension 
of phenomenally high temperatures and temperature gradients in the 
neighbourhood of the tool tip. Consequently, any theoretical 
treatment of the metal cutting process is bound to possess certain 
relaxations and simplifications. 

In the existing studies, simplifying assumptions are made 
with regard to the geometry of the deformation zone and the nature 
of stresses within it. These, along with the empirical 
determination of some of the process variables, have led to 
detailed calculations of the temperature field around the 
tool-tip. In the present study, the deformation mechanics has 
been treated in a fairly general manner and the temperature field 
predictions have been coupled with the deformation phenomena. Some 
of the involved features of metal-cutting, however, are not taken 
into account in the present study for the sake of simplicity. 
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2.1 PLASTIC DEFORMATION AND CHIP FORMATION IN ORTHOGONAL 
MACHINING 

In metal-cutting, the fact that the width of chip is large 
compared to its thickness, renders the problem two-dimensional. 
This is especially so, when the tool edge is orthogonal to the 
cutting direction. In the present work, a two dimensional 
orthogonal machining process under dry cutting conditions is 
considered as shown in Fig. 2.1. The tool is assumed to be sharp, 
with no built-up edge. The workpiece is taken to be ductile so 
that it produces a continuous chip. The tool and the portion of 
the workpiece material outside the domain shown in Fig. 2.1 are 
treated as perfectly rigid. Although the deformation of material 
in metal-cutting is similar to those of metal forming processes, a 
distinguishing feature is that the deformed material is separated 
from the parent workpiece in the form of a chip. Hence, the 
analysis of chip formation is tantamount to the study of plastic 
deformation phenomena in metal-cutting. 

The chip-formation can be described in terms of the following 
sequence of events. The uncut material in motion is first 
interrupted by the cutting tool. The tool is shaped in a suitable 
way to cause a large compressive stress on the material. After 
sufficient compression, when the resulting shear stresses in the 
material attain the yield value, plastic flow of material occurs, 
in the direction of the shear forces. Finally, the already 
deformed material particles on the rake face, are considerably 
slowed down due to the frictional stress prevailing between the 
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chip and the tool and also the work-hardening of the deformed 
particles. Further more, since extremely high-temperature 
conditions are produced near the tool-tip, the deformed material 
particles tend to stick partially to the rake-face. However, 
these particles also slip past the rake-face, due to the push 
generated by the subsequent elements of deformed particles. In a 
continuous operation, these sequence of events occur 
progressively to each successive element of uncut material. These 
events result in the continuous formation of the chip which is 
eventually separated from the workpiece. The plastic deformation 
and the frictional rubbing of the chip on the tool cause enormous 
amount of heat generation. The high temperatures resulting from 
such heat production, considerably influenced the deformation 
process and also affect the tool wear. It is therefore important 
to study both the mechanics and thermal phenomena during the chip 
formation process. 

Deformation mainly takes place in two regions, which are 
termed as primary and secondary zones. The region in which 
deformation occurs first is termed the Primary Shear Deformation 
Zone (PSDZ), wherein material is essentially sheared along a plane 
The PSDZ is denoted by the region ABCD in Fig. 2.2. The mean 
width of PSDZ is an important criterion which gives an idea of the 
rate of deformation and the amount of heat generated. The narrower 
the region, the greater the rate of deformation. Further, this is 
believed to be dependent upon the cutting conditions. For 
instance, the mean width of PSDZ is generally taken 
to be inversely proportional to the cutting velocity. 
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A portion of the chip which lies adjacent to the contact 
patch between the chip and the tool undergoes additional 
deformation because of friction at the rake-face (area PQRS as 
shown in Fig. 2.2). The name Secondary Shear Deformation Zone 
(SSDZ) is given to this region. Each layer in SSDZ moves at a 


di f f erent 

speed because of 

shearing in a 

direction 

parallel to 

rake-face . 

Although this region is small 

as compared 

to PSDZ, the 

amount of 

heat produced 

is believed 

to be of 

a similar 

magnitude . 

An interesting 

point to note 

is that 

the highest 


temperature occurs on the interface between SSDZ and the 
rake-face. The temperatures in SSDZ are higher than those 
securing in PSDZ mainly because the metal flow velocity is very 
low in the SSDZ and also the material which enters SSDZ has 
already been heated through PSDZ. At present, the thicknesses of 
SSDZ and PSDZ are estimated through empirical relations, which 
themselves are not on firm footing yet. It is definitely of 
interest to predict the extent of these zones, through a detailed 
analysis of the deformation phenomena. 

2.2 VI SCOPL ASTI CITY COMCEPTS 

It is well known that under tensile load, a ductile material 
deforms elastically upto a certain strain, beyond which it starts 
yielding. This limiting strain is known as the Elastic limit. 
Upto the elastic limit, the stress and strain in the material are 
related by the Hooke’s law. Beyond the elastic limit, when the 
material yields plastically, the stress (sometimes called as flow 
stress) becomes a function of the strain rate. Some of the 
properties of the solid undergoing plastic deformation are akin 
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to that of a Non-Newtonian fluid, in the sense that the 
relationship between the shear stress and the shear strain rate is 
non-linear. In fact, it is possible to define a pseudo-viscosity 
for a plastically yielding material which is the ratio between the 
shear stress and the shear strain rate. In the limit of a 
perfectly elastic or rigid solid, the viscosity assumes an 
infinitely large value. 


It is often assumed in the analysis of plastic deformation 
processes that the density of the material is unaltered during the 
process. The material incompressibility which is implied by such 
an assumption, is reasonable since yield occurs primarily due to 
shear. Under such conditions, material volume hardly undergoes any 
change. Therefore, the sum of all the diagonal components of 
strain tensor Ce^^), which represents the volumetric strain, is 
zero. That is. 


1 1 


11 


+ e 


22 


33 


= 0 


( 2 . 1 ) 


The strain tensor at a point itself can be expressed as 


e.. = -:re..<5.. + e.. 

ij 3 11 ij ij 


where 


e. . - strain tensor 

i J 

el . - deviatoric strain tensor 

J- J 

and 6..=1 ifi=j 

ij 

= 0 if i j 


( 2 . 2 ) 



15 


Assuminc material incompressibility to prevail, the various 
theories of plasticity developed so far attempt to relate the 
deviatoric part of strain (or its rate of change with respect to 
time) with the deviatoric stress. The deviatoric stress itself is 
defined in a manner similar to the deviatoric strain in eq. (2.2). 
Thus , 



o* . . - stress tensor 

I J 

O'. , - isotropic stress 
3 1 1 

O'! . ~ deviatoric stress. 

I I 

In addition to finding the relationship between the 
deviatoric stress and deviatoric strain, it is important todefine 
the criterion for material yield itself. An extensively ased yield 
criterion in the analysis of plastic deformation processes is due 
to Von Mises. According to this criterion the second invariant of 
deviatoric stress tensor (represented by J^) determines the 
condition of yielding. Mathematically, yield occurs when 


5 "ij 


( 2 . 4 ) 


vhere K is a material yield constant. 


Since material is said to be plastically deformed when 


yielding occurs, the total strain deviation is composed of two 
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parts i.e. elastic and plastic. Thus» one can write, 


e! . 
1 J 


i J i J 


(2.5) 


However, since the problems of metal forming and metal 
cutting involve huge plastic deformations, the elastic part of 
strain deviation is often neglected. This omission of elastic 
strain could introduce a maximum error of 2 to 3 percent which is 
permissibl e . 


Hence , 



( 2 . 6 ) 


of 
t . 


(p) 

The total plastic strain e: . occuring for a 

1 j 

time t is given by the area under the curve of 
Thus , 


loading period 

e. . versus time 
i J 


e 


(P) 

i j 


e^^^(o) + S' e . . dt 
o "■> 


(2.7) 


where e-^^o) is the plastic strain at the initial time t = o. 

1 J 

As a consequence of the principle of material incompressibi- 
lity (from eqns . 2.1 and 2.6), as 


o(P) - n 

e . , = o 


( 2 . 8 ) 


( P ) - 

the strain rate tensor e . in equation 

i J 

same as that of the pi astic deviation strain 


2.7 is therefore, 

tensor e.’ 

1 J 



1*7 


For varying loading 

. . ( P ) 

plastic strain e. . can 

i J 

the infinitesmal strain 
minute strain increments 


sequence, equation (2.8) for the total 
often be replaced by an algebraic sum of 
increments. The addition of all of the 
to obtain the total strain forms the 


essence of the Incremental theory of plasticity. 


2. 3 MATHEMATICAL FORMULATION 


A. CONSTITUTIVE RELATION FOR VISCO-PLASTIC FLOW 


In a metal cutting problem, it is conceivable to describe the 
metal flow, as that of a Von-Mises type of strain rate-sensitive, 
non-strain hardening, visco-plastic material, this, in turn, is 
akin to an incompressible (constant in density), non-newtonian 
fluid, as discussed earlier. An important aspect to be noted here 
is that the viscosity of such a "fluid” is dependent on the 
current local strain rates, total accumulated plastic strain 
undergone by the material and the local temperature. In the 
present analysis, however, work-hardening and thermal softening 
characteristics are assumed to be absent. Although the assumption 
has been invoked here for the sake of simplicity, a weak 
justification for the same can be offered based on the fact that 
work-hardening and thermal softening tend to cancel each other. 
Therefore, the material has been assumed to be purely 
rate-sensitive . For such a situation, the constitutive relation is 
expressed in the form : 


(2.9) 



is the 


derlvatorlc 


part of strain-rate tensor. The 


where e' j 


strain-rate tensor 

may 

itself be expressed 

in terms 

of 

the 

velocity gradients 

in 

the material. Its 

components 

In 

two 


dimensions (for orthogonal cutting) are given by 


du 


(2.10a) 
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(2.10b) 
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(2.10c) 


where u and v are the velocity components in x and y directions. 


The viscosity jj for a strain-rate sensitive material 

is taken to be a function of uniaxial yield stress of 

of the material (o- ) and the strain-rate invariant (e). 

y 

Mathematically, 




e) 


( 2 . 10 ) 


In the present study, the constitutive relations of the form 



(2 .11a) 


has been considered where y and n are the physical constants which 
define the visco-plastic characteristics of the material. They 
usually take values between 1 and 2LSee AppeiK'^ix b3 , 



Equation 2.11a can be looked at from another angle. As 
stated in Section 2.2, for a plastically deforming material, yield 
is a function of strain-rate (e) as well as the uniaxial yield 
stress (O' ) . This stress, often called as flow stress or the 
effective yield stress and denoted by is the term in the 
numerator of eq . 2.11a l.e. 


O* = O' 


y I /3 


l/n 


= fC^y , e) 


Rearranging terms in eq . 2.11a gives 


(2.11b) 


O' = -/3 fj e 


(2.11c) 


It may be noted that eq. 2.11c can also be derived from the 
Von-Mises criterion for yield and the stress-strain relationship 
of eq . 2.9. The strain-rate invariant e can be expanded in terms 
of the strain-rate tensor as 


In a 


e 



two-dimensional 


e 


a 


. . e, . 

1 j 1 j 

ituation. 


(2.12a) 

the above equation reduces to 
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(2.12b) 


In terms of velocity derivatives, e can be expressed in a form 
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(2 . 12c) 


B. GOVERNING EQUATIONS 

The necessity of considering the metal cutting problem in a 
different perspective from the earlier studies was highlighted in 
Chapter 1. Especially, the size and shape of the plastic zones and 
the temperature in the vicinity of the tool tip need to be 
predicted, without any adhoc assumptions or use of empirical data. 
Needless to say, this calls for a simultaneous application of 
mechanics of metal cutting and heat transfer theory. From the 
considerations of material balance, stress balance and heat 
balance, the governing equations in sbecfcdy state for the velocity, 
pressure and temperature fields are given by : 

(From Hass Balance) 

(2.13) 

7 . + f (From Stress Balance) 

~ ft; ~ 

(2.14) 

2 

and k 7 T + Q = p Cp (V .7 T) (From Energy Balance) 

(2.15) 


p( V .7 V ) = 



where, p = Density of workpiece material 
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= The Velocity Vector 
& = Stress tensor 
^ = Body force vector 

k = Thermal conductivity of work-material 
T = Temperature 

0 = Heat Generation Rate per unit volume 
Cp= Specific heat of workpiece material. 

FLOW EQUATIONS 


From equation (2. 4b), the stress tensor in the Momentum 
Equation can be expressed as 


a - -p I + d 



(2.16a) 


where, p is the pressure and, d is the deviatoric part of stress 

2b 

tensor . 


Also the deviatoric stress tensor d can be expressed in terms 
of the velocity gradients in the form : 


d = (V V + 7 V^) - I pi(7 .V) I (2.16b) 

Jy; 

T 

where V V represents the transpose of 7 V matrix. Simplifying 

/V rw 

eq.2.16b through the incorporation of material incompressibility 
condition and substituting for d in eq. 2.16a, one obtains 


■p I + (7 V 


T 

+ 7 V ) 


(2.16c) 
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Incorporating the above relation in the Momentum Equation 
(2.14), gives 


p(V . V V) 



+ /jCV V + V 


T 

V ) 


) 


(2.17) 


where the body force vector £ (which usually consists of the 
weight of material per unit volume) has been absorbed into the 
pr essure. 


COMPONENT FORM OF FLOW EQUATIONS 


As discussed earlier, the problem can be considered as a 
two-dimensional one, since the width of the chip is large in 
comparison to the thickness. The material balance equation for the 
two-dimensional situation can be represented as 


^u ^ £v 
Sx Sy 


O 


(2.18) 


where u, v are the velocity components in x, y directions 
respectively (See Fig. 2.3). 

The stress tensor o', in such a two-dimensional case, can be 
written in terms of its components as 


</ = ii + o- i] + cr ji + O' jj (2.19) 

.XX xy ■' yx ■' yy 

S|Sr 

Substituting for <y in the momentum equation (2.14), the x and y 

s» 

components of the momentum equations can be written in the form : 
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Y -MOMENTUM 
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(2 . 20b) 


Again the components <y , ty <y ^ 

XX' xy' yx' yy 

tensor, can be obtained from eq. (2.16c) as : 


the stress 
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(2.21a) 
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(2 . 21b) 


and O' 

yy 


-p + 



(2.21c) 


Substituting for these stress components in equations (2.20a) 

and (2.20b), the final forma of these momentum equations are given 

t>y 


P(u 
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(2.22a) 
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(-p + 2fj (2.22b) 
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ENERGY EQUATION 

In the energy equation, the important term to be carefully 
evaluated is the heat generation, since this can affect the 
temperature solution considerably. It is reasonable at this stage 
to postulate that the heat generated in metal cutting arises from 
two causes, namely : cl) plastic work in the PSD2 and (2) friction 
at the chip-tool Interface. The above assumptions are reasonably 
valid, because the^ are the most plausible ways in which mechanical 
energy is converted into heat. More importantly, it is not 
improper to assume that all of the work done is converted into 
heat; the two Important reasons for this are that the material 
density does not undergo a large change in value and also work 
hardening is not severe due to the high temperature conditions. 
Thus, the properties of the cut and the un-cut materials are not 
very different from each other. Keeping these factors in view, it 
appears reasonable to assume that all the mechanical power 
expended in machining is converted into heat. 

Since the plastic power expended is given by the product of 
the flow stress and the strain rate, the volumetric heat 
generation rate can mathematically be expressed as 

0 = T e 


(2.23a) 



zs 


where , 


'*■ - Effective yield shear stress or Flow stress 
and, e = Strain rate invariant. 

The Von-Mises criterion allows the shear stress to be 
evaluated in terms o£ the effective uniaxial yield stress of the 
material as 



=> from eq . 2.11c t = p e 


(2 . 23c) 


Substituting the above expression for t in eq . 2.23a, 
expression for the rate of heat generation per unit 
given by : 


the final 
volume is 


Q = ^ e C2.23d) 

VI 

The final dimensional form of energy equation can therefore be 
obtained from eq, 2,15 and 2.23d aa 


9 7 

ST ST 


VI 
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ax 
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(2.24) 



C. NON-DIMENSIOMALIZATION OF FLOW EQUATIONS 

As long as the number of parameters involved in a problem are 
few in number, it is realistic to stay with the dimensional form of 
variables, which enables ready comparison of pre.dLvcfcecl results 
with real values. However when the number of variables is large, 
analyzing the Influence of each parameter on the results becomes 
an arduous task. The subject of metal cutting is known to involve 
parameters, which exhibit very complex inter-dependences among 
themselves. To circumvent such hardships and facilitate a 
detailed parameteric study, it la desirable to non-dlmensionallze 
the variables of the problem with suitable scaling factors. 

A real difficulty which often arises during 
non-dimensional izat ion , is the proper choice of a scaling factor 
for particular quantity. Though some general guidelines do exist, 


it can be 

stated that the 

primary 

gui de 

in this 

process 

is 

intuition 

alone, since 

there 

is no 

unique 

way 

for 


non-dimensional izat ion . However, the beat possible way is to give 
the physics of the problem a critical examination and incorporate 
as much of the physical characteristics of the problem as possible 
while selecting the scaling factors. For the present problem, to 
begin with, the velocity components u and v are 
non-dlmensionalised with the cutting velocity (V); the distances 
with the depth of cut (t), and the pressure with the uniaxial 
yield stress (o* ). The scaling factors for the remaining 

y 

variables will be selected conveniently, during the course of 
non-dlmensionalizatlon itself. Thus, we define. 



(2.25) 




* * t -k 

X = x/t; y = y/t; V = v/V; u = u/V 


where the aeteriak denotes a dimensionless quantity. Using the 
above definition in material balance equation. 


dCu*V) d(v*V) 

+ = 0 ( 2 . 26 ) 
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one obtains 
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Similarly the x-momentum equation becomes 
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At this stage, the pressure and viscosity can be scaled as shown 
below : 
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(2 . 28b) 


Using the above expressions, the momentum equation can be 
simplified as 



Z2 



where o 
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(2 . 29b) 


Similarly, the 


y momentum equation can be 


reduced 


as 
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It Important to note that the Viaooaity e<,„atlon la alao to he 
normalized. Its non-dimensional isation gives the 
equation 


following 


fj = p/(o- t/V) = 


‘S' + f y— 1 

■/s t 
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Co-yt/V) 


(2 . 30) 
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NON-DIMENSIOMALIZATION OF ENERGY EQUATION. 


The energy equation in dimensional form as given by eq . 2.24 


is 


fh + 


+ ^ e 


^ , ar ^ ai . 


The normalization of temperature T in this equation can be 
done using some reference temperature » as 


ref 


(2.32a) 
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T f can bt conveniently chosen as shown below : 
ret ' 

Substitutinc for T and the other quantities in terms of their 
dimensionless counterparts, leadcto 


kT. 


ref 
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ax ay 
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- ref , 
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Multiplyine by t / (.pC ), the above equation reduces to 
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Taking o-y/(pCp) = and 


k/(pC^) = « (thermal diffusivity) J 


(2 . 32d) 


the final form of dimensionless heat balance equation becomes 
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where Vt/ot = Pe <Pec]et number) (E.33fa) 

2.4 BOUNDARY CONDITIONS 

The modelling of the boundary conditions is equally 
important as the modelling of differential equation themselves, 
since different boundary conditions may yield drastically 
different results. In fact, experience reveals that the toughest 
part of the problem formulation is the prescription of proper 
boundary conditions. Hence this section denotes due attention to 
the formulation of each boundary condition. 

To facilitate easy computation and to implement the boundary 
conditions in a systematic fashion, it is desirable to divide the 
boundary surface into individual segments, each with a different 
type of boundary condition. This division is essential 
particularly when the domain is complicated. For the present 
study, the total boundary has been divided into eight separate 
segments as shown in Fig. 2.3. 


The first three surfaces I, II and III which form the 


boundary in the workpiece material are essentially quite far 


off 




FIG. 2-3 PROBLEM REGIONS, SHOWING THE BOUNDARY CONDITIONS 
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from the actual recion of plastic flow and hence have prescribed 
velocity and temperature boundary conditions. The velocity 
components in the x and y directions are taken as unity and zero, 
respectively, while the temperature can be safely prescribed as 
the room- t emperatur e . 

The fourth surface, which is the already machined portion of 
the work piece may be assumed to be straight. Here the surface 
may be taken to be free of any boundary traction and the normal 
velocity may be taken to be zero. Thus, 

= 0 (2.34a) 

o 

and , 

V = o (2.34b) 

^n 

where F is the traction force and V is the normal velocity. 

^n 

Also, since the surface is in contact with the atmosphere, 
the convective boundary condition is applied for heat transfer in 
the form 

-k = h (T - T ) (2.34c) 

dn 

where n is the coordinate normal to the surface. 
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The fifth surface whose length is assumed as an input data 
(contact length between the chip and the tool), is one of the most 
difficult surfaces to be handled since a number of interesting 
events occur on this surface. Since the chip is in perfect contact 
with the tool over this surface, the normal velocity on this 
surface can be prescribed to be zero. Thus, 


= 0 (2.35a) 

A point to be hifihlifihted for this surface is that unlike all 
the other sefiinents, this experiences boundary tractions arising 
from the frictional and compressive stresses due to direct contact 
with the tool rake-face. However, the normal traction need not be 
specified since has been prescribed as zero. As regards 
frictional stress, it is assumed that this stress is given by the 
product of the local yield stress and a constant friction factor. 
Thus, 


T ^ = m 


■/3 


(2.35b) 


= > from eq . 2.11c = m e 


(2 . 35c) 


Assumption of a 
contact length is 
However, for want of 


constant friction factor for the entire 
questionable, as discussed subsequently, 
a better estimate of variation of m along the 


contact patch, a constant value has been assumed. 



2 ) 5 “ 


Re^ardinfi the heat transfer conditions, the heat generated by 
friction at the surface has to be accounted for. This is a very 
important condition, as it leads to the occurrence of a high 

temperature region over the contact patch. In the present study, 
it is assumed that the heat generated due to frictional rubbing is 
shared between the chip and the tool in the ratio of the 
respective thermal conductivities. Thus, 

k 

-“SK " e; (2.35d) 

where k and are the thermal conductivities of the chip and the 
tool and is the rate of heat aeneration per unit area due to 

friction. From basic mechanics the frictionl power dissipation is 

given by the product of shear stress and velocity. Ilathematically , 

Qf = Tj (2.35e) 

where t, = Effective frictional shear stress (=ro— or «uje ) 

' 1/3 

V.J. = Tangential velocity of slip 

0^ = Rate of Frictional Heat Generation. 


The tangential slip velocity V.j, can, in turn, be evaluated as : 
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= V cosoi + u sinot 


where a = Rake Anfile of Tool. 


(2.35f ) 


Conibinin£ equations 2.29c, 2.29d, and 2.29e the final form of 
heat transfer equation for the fifth surface is : 


_ k Jnky 
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(2 .35g) 


For the sixth and eight surfaces, the conditions of zero 
normal velocity, zero surface traction and convective heat loss to 
the atmosphere are applicable. Therefore, on these surfaces, 


V = 
^n 


0 
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F 
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(2 . 36b) 


and , 
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(2 . 36c) 


For the seventh surface, which is on the chip exit side, the 
following conditions are applicable. From a mass balance between 
the uncut and the cut material. the chip velocity can be 
estimated. Since the chip material is assumed to become rigid once 
again far away from the tool tip, the velocity can be taken to be 
equal to V on this boundary. The velocity components are then 


given by 



57 


u = V sino< 
c 


(2.37a) 


V = V cosot 
c 


(2 . 37b) 


where 


V = V i— 
c t 

c 


(2 . 37c) 


For heat transfer, there is no suitable boundary condition. 
It is assumed that the temperature of the chip becomes invariant 
in the chip-flow direction. Thus, 


Sn 


= 0 


(2 . 37d) 


NON-DI KENSa; ONXLI ZATI ON 

The non-dimensional form of the boundary conditions are: 
For convective heat-loss on boundaries IV, VI and VIII, 


= > 


dT 

Sn 
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h (T - T^) 


where h = 


h 

WtT 


and 


T, = 


ref 


(2 . 38a) 


(2.38b) 


On the chip-tool interface (surface V) 


* * 


fe T V 

Sn t 


(2.39) 
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On surfaces I, II and III, 


* * 
T = 


(2.40) 


At the far end of the chip (surface VII), 


Sn 


= o 


(2.41) 


The dimensionless velocity conditions are : 


★ ic 

u = 1, V = o on surfaces I, II and III (2.42) 


= o for surfaces IV, V, VI and VIII (2.43) 


The chip velocity V^ can also be non-dimensional ized as 


* 

V = 


(2 . 44a) 


resultin£ in 


* t 

u = V cosoi 
c 

* * 

V = V sinoi 
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y For surface VI I 


C2.44b) 


The frictional stress condition on surface V becomes : 


* 1 + B* (e )^^" 

■"f = 




(2.45) 
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The dimensionless ftovernin£ equations and boundary conditions 
formulated above have been solved by the application of the Finite 
Element Hethod. The details of the solution procedure are 
described in the next chapter. 



CHAPTER 3 


FINITE ELEMENT ANALYSIS 

In many engineering problems, it is not always possible to 
obtain a closed form exact solution. Naturally, recourse needs 
to be taken towards approximate solutions in such cases. 

Fortunately, with the advent of digital computers, the effective 
application of approximate numerical techniques in almost all 
branches of engineering has widened the scope of software solu-ttoris 
by leaps and bounds. The finite Element Method adopted for the 
present study is one such powerful technique for obtaining 
numerical solutions to difficult problems. 

3.1 WEIGHTED RESIDUAL METHOD 

The differential equations corresponding to an engineering 
problem are required to be satisfied at every point in the 
solution domain. However the numerical solution for the problems 
may not satisfy the governing equations exactly. Instead, it may 
satisfy the equations approximately, leaving a small non-zero 
residue at moat of the locations in the solution domain. 

* 

Denoting the exact and approximate solutions by * and S , the 
substitution of these solutions into the governing differential 
equations yields 

LC*) = 0 Exact Solution 

= R Approx. Solution 




LC») 


0 


Exact Solution 


(3.1) 


LcaT) = R 


Approx. Solution 


where L is the differential operator of the probleoi 


From the above expressions, it is obvious that the smaller 
the value of residue, the closer the approxim?ite solution is to the 
exact one. An Alternate way of looking at it is that the larger 
number of points at which the residue is small, the closer are the 
solutions 4 and # . Since, limitations exist on the computational 
time and effort, it may not be possible to increase this number 
beyond a certain value. Therefore, methods have been developed by 
which the residue la minimized in an average (integrated) sense 
over the solution domain. In order to pin-down the residue to 
small values at chosen number of points in the domain, approximate 
weighting functions are multiplied during the minimization of the 
residue In the whole solution domain. 


The above discussed objectives are achieved 


J R dv = 0 

D 


for 1 = l,..n (selected nodes) 


where are chosen weighting functions and D 


by setting, 

(3.2a) 

is the solution 


domain. 
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Substituting eq . 3.1 into eq. 3.2a, it gives 


J 


u 


L(* ) dv = 0 

£of 1 = 1 , . .n 


(3.2b) 


Equations (3.2b) form the basis of all Weighted Residue 

methods, for solving differential equations. The manner in which 

the weighting functions are selected and the approximate solution 
* 

5 is defined D, leads to many different weighted residual 
methods such as Ritz Method, Galerkin Method, Collocation Method 
and the method of Least Squares. The Garlerkin’s method has a 
general applicability and hence been adopted here for the Finite 
Element solution procedure. 


3.2 FINITE ELEMENT PROCEDURE 


The finite element procedure as the name indicates, implies 
that the solution domain is divided into many finite sub-domains 
(or elements). This discretization is advantageous particularly 
for complex domains since the focus of attention now is reduced to 
that of an element of chosen shape. Considering the element as a 
building block, complicated domain shapes can be represented, even 
by the use of elements of a particular shape. Thus the analysis 
is very much generalized and computation becomes easier, natural 
question which then arises is how these elements are related or 
clubbed together. The answer lies in the fact that since the 



4 ?) 


va.ria.ble8 vary continuoualy acroas the eletoenta, a group of 
eleoienta could be aaaembled to form a larger domain. 

The ae(]uence of procedurea in any FEM aolution atrategy could 
be outlined aa followa. 

1. Diacretizat ion of domain into amall elementa. 

2. Derivation of elemental propertiea. 

3. Grouping the elementa into a global aaaembly. 

4. Incorporation of Boundary Conditlona. 

5. Solving the reaultlng global Matrix Equationa to obtain the 
field varlablea. 

6. Any proceaalng that can further be done. 

The above mentioned atepa are discuaaed in detail later. The 
advantage in aplitting the domain into many elementa ia that the 
field varlablea can be Interpolated using standard interpolation 
functions within each element. For this purpose, nodes are 
selected within each element and the unknown variable ia expressed 
in terms of the values of field variable at the selected nodes. 
The interpolating functions for an element are also called as 
Shape Functions. Every node has a shape function whose value la 
unity at that particular node and zero at all the other nodes of 
that element. This la a common property of all the shape 
functions, although in their form they may differ from each other. 


Mathematically, the value of the field variable within an 
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element can be expreaaed as 


* 


n 



*1 


coKcve ate nodni valu.GS (3,3) 


Since such representations are possible for all the elements in 
the domain, it can be said that the field variable la known in the 
entire domain. Using standardized shape functions, elemental 
shapes and nodal locations, the elemental properties are easily 
described. These, In turn, aid the computation of the solution 
on a computer. 


3.3 TYPES OF FLOW FORMULATION 


The nature of a flow problem and the flow quantltltes of 
Interest (say velocities, pressure, vorticity etc.) Influences the 
solution approach. Based on these considerations different types 
of flow formulations have been developed which are 
explalTfied brl ef ly here. 


The Stream-Funcblon approach, first proposed by OSLON 
reduces the number of variables in the flow problem to a single 
variable. This new variable stream function 4' is defined in such 
a way that it automatically satisfies the mass balance equation. 
Also the Mbmentuffl Equations are reduced In terms of the stream 
function. But this process results in a higher-order differential 
equation (fourth order) and demands higher-order Interpolations 
within the elements which Is disadvantageous in some cases. More 
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importantly, to obtain the velocities results have to be further 
processed. In the present problem, since velocities are directly 
used in the energy equation, this approach is not adopted. 

In the Stroam-FuncLlon Vorblciby approach an additional 
variable called vorticity (represented as is defined and the 
three equations of mass and momentum balances are reduced to two 
2nd order equations in 4' and Ci>. As mentioned above, the 
primary interest in the present problem being velocities, this 
approach is also not suitable. 

The third method called the Velociby-Pressure involves the 
direct solution of governing equations. There are many advantages 
for this approach such as applicability to three-dimensional 
flow, easy incorporation of pressure and velocity boundary 
conditions and so on. Also zeroth order continuity being 
sufficient for the elemental interpolation functions, this 
approach requires less computational time than the other 
approaches. Often, to eliminate some of the numerical difficulties, 
the pressure variable is elimlnted by taking mass balance as a 
constraint and Introduc into the momentum equation. This is 
done by what la called as the penality function procedure. Bat this 

needs that the pressure boundary conditions be modified 
accordingly and hence is believed to be unsuitable for the present 


problem. 
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3.4. APPLICATION OF FEM 


A. DOMAIN DISCRETIZATION 


Discretization of the domain involves division of the 

solution domain into small elements. For plane two-dimensional 
problems, these are n-sided polygons and hence many types of 
elements are possible, like triangular , quadrl lateral etc. The 
domain discretization, seemingly simple, can affect the accuracy 
of the solution considerably, if done Improperly. This is 
particularly so when the domain is of a complex shape. 
Isoparametric elements (so-called because the field variable and 
the spatial coordinates can be defined with the help of the same 
interpolation f unctions Joffer many advantages and hence employed 
in the present study (shown in Fig. 3.1). Elght-noded linear 
quadrilateral elements have been used, as shown in Fig. 3.2, in 
creating the finite element mesh. 


B. DERIVATION OF FINITE ELEMENT EQUATIONS 

The finite element equations for plastic flow of metal, are 
better derived from the basic vector equation (eq. 2.13) so that 
the handling of the traction boundary conditions at the chip-tool 
interface could be easily implemented. 

FLOU EQUATIONS 

The vector Momentum Equation eq. 2.14 is given by 



Eight Noded Element Showing the Local 
Co-ordinate System. 
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p(V.VV) ~ ^ . Cf - f = [0] 
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By Galerkin technique, the vector residue 
two-dimensional cutting can be formed as 


equation for 


II I pV. 7 V - 7.0- ~ ^ I 


dx dy = [0] 


(3.4a) 


= > 


II (pV . 7 V)dxdy - || (7.<y) dxdy - || f dxdy = [0] 


The term f in eq. 3.10 can be 
redefining pressure as 
V P’ = 7P + f 


(3.4b) 

included in the pressure by 

(3.5) 


The term JjNl (7.0-) dxdy can be written by the Ueak 
Formulation as 


JJ (7.0) dxdy = II 


D 


7. (N^ 0-) - 7 .o) 


dxdy 


(3.6a) 

Expanding the first term in parenthesis by divergence theorem, the 
above equation 3.6a reduces to 


11 ' 11 7 . o dxdy 


(3.6b) 
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where In the first term gives the contributions due to the 
specified boundary conditions. The original momentum residue 
equation 3.4b can now be rewritten as 



p (V. 7) V dxdy + JJv Nj 


= # n . c dl 


a 


dxdy 


(3.7) 


Expanding ct into its components, the LHS term 7N. .«>■ can be written 
as 


7N. .O' 
1 ^ 




j 


X-Momentum 


(3.8) 


Putting the above term of eq . 3.8 into the momentum equation 3.7, 
the X and y component equations can be derived as 



dxdy 


= # N 


i 


(n .O', i ) dl 



(3.9a) 


and 
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y-Moroentum 


JJ*”' [■ : 

If P [' 


dv'j 

dyj 

dxdy 


1 


r 

Sv 

dN. 



- p * 



k. J 


O', j) dl. 



dxdy 


(3. 9b) 


In the above equations, the inertial terms have n 
quasi-linearized by using the velocities corresponding to prev a 
iteration (or initial guess) values during the iterative solution 
of the problem. ' For this reason, these quantities have been 
denoted as u and v. 

The field variables have been interpolated within each 
element as shown below 



(3.10a) 


Similarly, the Spatial Coordinates of 
isoparametric element are defined as 


the boundaries of 


the 
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X 


n 



n 

y = ^ Nj yj 

1 =1 


(3.10b) 


where n = 1,8 
01 = 1,4 


It may be noted that pressure is defined only at corner nodes 
of the element and hence needs only a linear interpolating 
function . 


All the other variables (u, v and T) have been interpolated 
using quadratic shape functions. The lower order interpolation 
for pressure has been found necessary to enhance numerical 
stability [10]. 


For easy computation, the interpolation function are defined in 
terms of local coordinates as 
Corner Nodes 


Ni = 4 

Mid-side Nodes 


(1 + T7 - 1) (i+Vin) (3.11a: 


Ni = 2 ~ ^ ^ 

1 2 
Ni = (1 + ? ^ ) (1 - 


(3.11b) 


107.912 



Substituting the expressions for the interpolation of u, v 
and p Into the Momentum equatlona 3.9a and 3.9b givea the final 
non~dimena ional form of x and y momentum equatlona aa 


X-Momentuffl 
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(3.12b) 



The welghtlne fuftctlona for the mass balance equation have, 
been choaen to be the linear ahape functlona uaed to interpolate 
preaaure within an elenient. Thla ia done in order to get a 
preaaure field that ia coaiparable with the incompreaslbillty 
condition. Subatitutlng the interpolation expreaalon for the 
velocity variable (eq. 3.10a), the final form of the maaa balance 
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equation is 


n 




J 




[Uj] 


11 


dN 


— 3x 


[Vj] = 0 


C3.15) 


ENERGY EQUATION 


The weighted realdue form of the energy equation (eq. 2.24) 


ia 





2 '' 


+ 

9x dy 


dxdy + 



Be e dxdy 
y3 




dxdy 


0 


(3.16a) 


Since by Green'a theorem 



SI 

* “i 


dl 



(3.16b) 


ff 

ST 


-1 

JJ 1 *^* 

Sx 


dy J 


dxdy 

expanding the temperatlre within an element as T 
the above equation reduces to 


E »j Tj. 


rr 1 f !!!i !!Lt * “j 1 

JJ I ^ ^y J 


r - 1 

N. Pe N.U. ^ + N^V. ^ 

I [^kk^x kk^yj 


dxdy 


[Tj] 


= II 


N. Pe — e dxdy + # N. 4— 
i yj L Sn 


(3.17) 


where the first term on RHS la the heat generation contribution 
and the second term la from the heat transfer boundary conditions. 

The weighted residue equations for momentum, mass and energy 
balances, lead to matrix equations in terms of the nodal variables 
of velocity, pressure and temperature. The contributions to the 
matrix equations from each element can be evaluated by applying 
the momentum and energy equations at all eight-nodes of an element 
and the mass balance equations only at the four corner nodes of 
the element. Thus, for each element, the coefficient matrix 
contribution la a (28x28) matrix, while that for the right hand 
aide la a vector of size (28x1). 
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After computation of the element matrlcea, 
conditions are incorporated which la explained 
section. Later the elemental matrices and vectors 
into global matrices, vectors and the resulting 
solved iteratively. 


the boundary 
in the next 
are assembled 
equations are 


C. BOUNDARY CONDITIONS 


For the present problem, the domain as well as the boundary 
conditions are complex. To facilitate numerical computation, 
certain assumptions had to be made with regard to the shape of the 
chip f ree-surface , the contact length, the friction factor on the 
rake-face etc. In this section, the implementation of the 
boundary conditions on all the eight boundary surfaces through the 
evaluation of the appropriate boundary Integrals or by the 
prescription of the known values of velocity, pressure etc. is 
explained . 


Surfaces I, II and III have prescribed boundary conditions 
for the variables u, v and T. These conditions are Implemented by 
incorporating them as nodal equations for the corresponding 
variables, in the final matrix equations. 

The normal velocity (V^ = 0) la applied on surfaces IV, V, VI 
A VIII, since there is no flow in a direction normal to these 
This condition la satisfied by expressing the normal 


surfaces . 



componenta u 


and V and the 


velocity in terms of the velocity 
direction coainea o£ the normal at 
velocity at the auriace requirea 


that node. 


Thua, normal 


= u Coa © i - V Sin © J = 0 (3.18a) 

where & la the angle made by the normal with the x-axla and 1, J 
are unit vectora in x and y directlona respectively. 


Again, expanding the tecma u and v aa in eq. 3.10a, In terma 
of the nodal velocity values and interpolating functions., the 
above equation can be rewritten as 


V 

^n 



u 


Cos & 


i 



Vj Sin e j 


0 (3.18b) 


where u^, v^ are the nodal velocities of the element which lies on 
the concerned boundary. 

The angle made by the surface normal with the x-directlon can 
be given as input or can be evaluated from the element shape 
(specified by nodal coordinates) using the iso-parametric nature 
of the element. 

On surface VII, which is the far end of the chip, the 
velocity is equal to chip velocity V^. The velocity components u 
and V on this boundary are given by 
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u = V Sin ot 
c 

V = V Cos ot 
c 

where <x = rake angle of the tool. 


(3.19) 


On surface V, In addition to zero normal velocity, the 
traction condition due to Interface friction Is to be Implemented. 
As shown by equations 3.12a and 3.12b, the traction contributions 
to the X and y momentum equations are respectively 


X-Momentum 


J Nj dl and 

V 


y-Momentum 


I "i S 

V 


(3.20a) 


(3.20b) 


where t^ and 
directions . 
given by eq. 
the current 


t are tractions per unit area, in the x and y 
On surface V, the shear stress can be evaluated (as 
2.35c in terms of the friction factor, viscosity and 
strain-rate. Thus, in terms of t^, t^ and t^ reduce to 


the form 



X 


( 3.21a) 


5 3 


t 


T ^ Cin <x 


and ty “ f " 


(3.21b) 


Surfaces IV, VI and VIII are free as no traction is applied 
on any of these. Therefore, the traction integrals for these 
boundaries are identically zero and only the zero normal velocity 
condition is applied on these surfaces. 


As regards the heat transfer boundary conditions, these 


Involve either given 

temperature 

or 

flux or 

convective 

conditions. Surfaces 

I, II and III 

have 

prescribed 

temperature 

whose implementation 

has already 

been 

described . 

On other 

surfaces, the relevant 

boundary heat 

flux 

Integral term 

P ai 

J ''i ^ 

dl 
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a 


needs to be evaluated. 

For thesur faces IV, VI and VIII which have convective heat 
transfer b.c. , the overall heat-transfer coefficient is to be 
evaluated from forced convection correlations for air 
present study a simplified expression of the fotm 

so-/~w 


h 


(3.23) 



GO 


h&a been chosen, where V la the cutting velocity in m/a 
The boundary integral term of eq . 3.22 

f dT 

J “i ^ 

O 

can be rewritten using the dimenaionlesa boundary 
given in 2.38a in a form 


J" 


dT 
i in 


dl 


- Jn, 


h (T - T^) dl 


On the surfaces IV, VI and VIII expanding T from eq . 3 


= I 

J = 1 


the expression for convective heat loss becomes 




h (T - T^) dl = 


|(N^ h Nj 


) dl 


tT^] 


JCN. h T, 


)dl 




In the above equation, the first term on right 
to the left hand side of the overall matrix equation 
second one contributes to the right hand side vector. 


condition 

C3. 23a) 

10a as 


(3.23b) 


contributes 
while the 
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On the chip-tool interface (surface V), the condition of 
frictional heat generation haa been taken into account. As 
described in Chapter 2, the heat flux entering the chip-side is 
calculated through the boundary Integral 

r aT r 

j “i ^ J “i Vydl (3.24a) 

V V 

where the dimensionless shear stress itself is evaluated from 
eq . 2.45 and the slip velocity la given by 

= V Cos oi + u Sin a (3.24b) 

The frictional heat flux contribution given above is added to the 
right hand aide of the nodal equations of the temperature 
variables, on surface V. 

D. ELEMENT ASSEMBLY 


In the previous section, the elemental contributions to the 
left hand side coefficient matrix and the right hand aide vector 
were discussed in detail for the eight-noded lao-parametr ic 
quadrilateral elements. These elemental contributions are 
assembled into a global matrix equation by adding the entries 
corresponding each nodal variable appropriately. 



The boundary Integral contrlbutlne on both left hand and 
right hand aides parts of the matrix equation are also 
incorporated as discussed in the previous section. This assembly 
procedure results in a matrix equation of the form 

[A] [XJ = [B] (3.25) 

where [A] = Coefficient matrix 

(X] = Solution vector 

[B] = Right hand side vector. 

The problem at this stage is set for solving the above matrix 
equation by any standard technique available. 

E. MATRIX SOLUTION TECHNIQUE 

The selection of the matrix solution technique depends upon 
the matrix characteristics such as diagonal dominance TRID 
structure or pentadlagonal structure etc. 

For the present problem an iterative technique is used for 
matrix solution based on the Frontal method. The Frontal method 
offers many advantages such as small use memory space and high 
speed of computation. It utilizes the fact that when all the 
contributions for a particular node and for a particular variable 
are over during element assembly, it could be transferred to the 



disk inernory. 


On this basis it retains in its ineinory only those 
variables which are yet to be assembled. After all the elements 
are assembled, the variables are solved by Gaussian elimination. 
Since the governing equations are highly non-linear an iteration 
procedure with successive under-relaxation has been employed. 


3.5 PROGRAMMING 

Computer Programming has been one of the challenging tasks of 
the present work. The challenges arose both due to the complex 
nature of the problem and the boundary conditions. A 
brief description of the programming aspects are discussed here 
for the benefit of interested readers. 

Appendix C gives the flow chart of the program, emphasising 
mainly what each of the subroutines do. The Main Program calls 
the Subroutines of DIMENS^ DINPUT, ITERAT which in turn call the 
other subrout inesCAfP®''^’^* • 

Subroutine DIMEMS returns the values of the array dimensions 
which remain unchanged through out the program. Subroutine DINPUT 
reads and returns all of the required input parameters such as the 
axi -symmet r i c flow indicator, the number of initial and boundary 
conditions, body forces, tolerance, relaxation factors, coordinate 
information, scale factors, element connectivity data and the 
prescribed values of parameters like the Ambient temperatures. 
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material properties etc. The whole input data for the problem is 
non” d i m ens i ona 1 i z ed within the program. Also some of the fiiven 
input is checked for correction by two routines named DIAGNl and 
DIAGN2 . 

Subroutine ITERAT calls the FRONTS routine for solving the 
assembled matrix equation. This routine begins the problem with 
guessed initial values for all the nodal variables and relaxes 
them for the next iteration in the subroutine TOLREL if the 
solution is found not to have converged within the prescribed 
tolerance. This routine also calls the WRITER subroutine Ci) to 
write the nodal values of the variables if the solution has 
converged and (ii) to write the largest relative change in the 
variables that has occured during a particular iteration, if the 
solution has not converged. 


The FRONTS does 

the 

assembling 

of 

a 

group 

of element 

and 

finally solves for 

the 

variabl es 

by 

Gaussian 

elimination 

and 

back- substitution. 

It stops the execution 

of the 

program 

under 


two conditions. Firstly, when the front width prescribed is found 
to be too sail or if a pivot value is foud to be smaller than the 
set value (usually ste 10 For assembling it calls the MATRIX 
subroutine which calculates the element fluid matrix FLUMX. This 
fluid matrix consists of the elemental contributions of the 
governing equations to form the coefficient matrix. Similarly, 
the elemental right hand side vector is also calculated in the 
subroutine MATRIX for each iteration, the subroutine DRIVES is 
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called by MATRIX, which in turn calls the subroutines SHAPES, 
DJACOB & SHAPEA to calculate the shape function and their 
derivatives . 

The boundary conditions are implemented in the BOUCON which 
is called by MATRIX routine before it returns to PROMTS routine. 
This is called by identifying the surface, element and the side 
which are forming the boundary. The relevant boundary conditions 
from among the zero normal velocity, convective heat loss, the 
frictional heat flux and the frictional shear stress are 
identified appropriately and the elemental contributions are 
evaluated. The contributions to right hand side vector are 
calculated at local nodes and are assembled in the global vector 
both in MATRIX and BOUCON subroutines. 



CHAPTER 4 


RESULTS AND DISCUSSIONS 

The velocity and the temperature fields in the vicinity of 
the tool edge during metal-cutting have been obtained by the 
finite element solution of coupled stress balance and heat balance 
equations. The vork material has been taken to be purely 
strain-rate sensitive and obeying the Von-Hlses yield criterion. A 
parametric study has been conducted to highlight the effects of 
some of the Important process parameters such as Cutting Velocity, 
Depth of Cut and Chip Thickness Ratio. The properly values and 
the range of parameters used in the study are summarized in 
Appendix A. 

From the velocity and temperature fields some of the other 
features of interest such as the shapes of primary and the 
secondary deformation zones, the location of the shear plane and 
the mean width of PSDZ have been obtained. The effects of cutting 
parameters upon the temperature variation over the chip-tool 
Interface region, have been studied. This la very important since 
the highest temperatures occur in this region, causing a large 
tool wear. The structure of the primary and secondary deformation 
zones have been further analyzed by plotting the strain-rate 
variation In these zones. Comparisons with earlier works have 
been provided where available. 



4.1 DEFORMATION PATTERNS AND ISOTHERMS 


Figures 4.1 to 4.6 give the Deformation patterns and the 
corresponding Isothermal patterns for six different sets of 
cutting conditions. The data for the cutting parameters have been 
taken from [16]. In order to obtain the deformation patterns, the 
following procedure has been employed. The strain-rates (e)at 
different points were first normalized with respect to the maximum 
strain rate occurlng in the solution domain. Points having the 
same value of normalized strain-rate were connected by smooth 
curves and the curves of 3%, 4% and 5% strain-rate (as compared to 
the maximum) were thus plotted. These constant strain-rate curves 
provide a clear picture of the primary and secondary deformation 
zones (PSD2 and SSDZ). The isotherms have been plotted by 
Interpolating the available nodal temperatures. 

In Fig. 1(a), the deformation pattern la depicted with the 
help of 1%, 2i, 34, 44 and 54 strain-rate curves. The above 
values were used to obtain an idea of the approximate cut-off 
value that can be taken to demarcate the primary shea/ deformation 
zone. It can be seen that curves of 34, 44 and 54 cut-off values 
are relevant in this context and hence been used for other figures 
(Fig. 4.2 to Fig. 4.6). Indeed it appears that the 54 curve 
possesses a good correspondence with experimentally observed 
shapes, although an appropriate cut-off value can be proposed 
only after a detailed comparison of the computed results and the 
phot crographa £or alfnllar cutting condltlona. 
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The followlnfi interesting features can be noted from 
Fig. 1(a). Firstly the constant strain-rate curves on the uncut 
material side are more closely placed than the cut material side 
(in the chip). The reason for this appears to be that theincomlng 
material undergoes deformation due to the shearing action in the 
vicinity of the shear plane, caused by the large shear stresses 
prevailing there. Uithln the chip, however, the deformation is 


due to the curling 

of the chip and the 

compressive 

force 

of 

the 

tool on the chip. 

Since one of the 

surfaces of 

the 

chip 

is 

unrestrained, the 

deformation in the 

chip occurs 

over 

a wider 


region. It is to be noted, however, that the deformation near the 
free surface does have some dependence on the prescribed shape of 
this surface. The figure shows that a small region in the 
workpiece directly below the tool-tip is also getting deformed, 
because this la the place where the material is getting separated 
from the workpiece. Commensurate with the high stress intensity 
here, the strain-rate curves are also closely spaced. The 
secondary deformation zone (SSDZ) is observed to be approximately 
triangular in shape which is usually seen in experiments too. As 
expected this zone extends only till the contact length for the 
chip with the tool. Close spacing of the strain-rate curves is 
seen for SSDZ also, as discussed in the other zones, due to large 

stresses . 

For the temperature field also (in Fig. 4.1(b)), the trends 

to bo roollotlc. Tho «axlmo<o tomporaturo occata aoaawhara In 
the region where the chip la In contact with the tool. Thla haa 



well and is explained in 


t>een noted by earlier researchers as 
terms of the fact that in the SSDZ further rise in temperature 
occurs due to frictional heatlnfi beyond the level of heating in 
PSDZ. For the same reason, the gradients in the secondary 
deformation region are the highest. The lowest temperatures in 
the chip occur in the middle rather than at the free surface in 
accordance with the deformation pattern. Bulk of the material is 
at room temperature not being influenced much by the heat 
generation. This la because of the large amount of heat removal 
by the chip. The temperatures in the primary zone are far less 
when compared to those occuring in the secondary zone. This also 
establishes that the frictional heat generation on the chip-tool 
interface is very large as compared to the heat production in 
PSDZ. The workpiece material below the tool-tip is also getting 
heated because of high deformation in this region. The occurence 
of very high temperature on the rake face of the tool has an 
important bearing upon the tool wear. The coupled analysis 
presented here can therefore be utilized to calculate the stresses 
and the temperature on the tool surface which, in turn, can be 
linked to tool wear. 

Figures 4.1 to 4.6 also show the affects of the cutting 
conditions on the deformation and isothermal patterns. 
Particularly, Figs. 4.7 to 4.9 respectively show the effects of 
cutting velocity, chip thickness ratio and depth of cut on the 
temperatures at the chip-tool interface. The deformation regions 
(both primary and secondary) are wider for larger cutting 




FIG.4-i(a) DEFORMATION PATTERN 
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FIG 4-E(a) DEFORMATION PATTERN 
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FIG.4-3(b) ISOTHERMAL PATTERN 
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FIG.4.6(b)IS0THERMAL PATTERN 
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Rg. 4.7 Effect of Cutting Velocity on the 

Temperatures at the tool-chip interfoce. 
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Rg. 4.9 Effect of Depth of Cut on the 

Temperatures at the chip— tool Interface. 
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velocity. Conaequently the temperatures are also higher, 
especially at the chip tool Interface also seen from Fig. 4.7. It 
can be seen from Fig. 4.8 also that the chip thickness ratio has 
only marginal effects on both the temperature fields and 
deformation patterns. But the Isothermal patterns are very much 
influenced by the depth of cut as also seen in Fig. 4.9. However, 
the primary deformation zone temperatures are affected much more 
than the secondary zone temperatures. This could be because the 
deformation may not occur uniformly across the shear plane, when 
the depth of cut is Inci^eased. The Isotherms also seem to take 
Irregular shapes in the case of large depth of cut, due to the 
same reason. 

4.2 IDENTIFICATION OF THE SHEAR PLANE AND VARIATION OF PROPERTIES 

ALONG AND ACROSS THE PLANE 

The shear plane could be identified using the available 
strain-rate distribution. The constant strain-rate curves for 
higher and higher values were plotted on either side until they 
approached very close to each other becoming almost parallel 
lines. This . « shown in Fig. 4.10 where the curve IV has a 
strain-rate (e) of 10% of the maximum. The plane AB shown in the 
figure can thus be said to represent the shear plane. 
Interestingly, the shear angle obtained in this fashion matches 
well with the measured one (given in table). 
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Figure 4.11 shows the same case with the shear plane and the 
primary deformation zone. A cut-off value of 5% has been used to 
Indicate the boundaries of PSDZ (Curve III). To obtain the mean 
width of the PSDZ, normals were drawn from the shear plane to 
Curve III and the average length of these normals was taken as the 
mean width. 

Often, the strain-rate and temperature variations across and 
along with shear plane are of particular Interest. In Figs. 4.12 
and 4.13, the strain-rate variation across and along the shear 
plane have been drawn. As can be expected, across the shear plane 
the strain-rate rises rapidly, reaches a maximum somewhere near 
the shear plane and falls rapidly again. Also the strain-rate is 
seen to be very high near the tool tip and small in the middle of 
the chip. It increases slightly as compared to the central 
portion, near the free surface. However, the range of variation in 
the strain-rate is leas along the plane when compared to that 
across the shear plane. The same trend is observed for the 
temperature variation shown in Figs. 4.14 and 4.15. 

4.3 VARIATION OF STRAIN-RATE NORMAL TO THE CHIP-TOOL INTERFACE 

In Fig. 4.16 and 4.17, the variation of strain-rate and 
temperatures in a direction normal to the rake face are shown for 
three locations, viz. (i) near the tool-tip, (ii) middle of 
contact length and (iii) end of the contact length. As can be 
seen the strain-rate and the temperature decrease rapidly after 




Fig, 4.16 Venation of Strain Rate Normal 
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the aecondary deformation zone. An idea of the structure of SSDZ 
can be obtained from these figures also. Similar to the 
determination of the width of PSD2, the dimensions of SSDZ can 
also be obtained by drawing suitable cut-off curves beyond which 
the strain-rate and temperature gradient become very small. 

The results of the present study confirmed all the well-known 
trends in metal cutting. However, one peculiar feature of the 
prediction obtained in this study is that local maxima in 
temperature are seen at two different locations in the chip-tool 
Interface which has not been observed by the other researchers. 
One of these maximum values occurs near the tool-tip while the 
other is close to the location where the chip separates from the 
tool surface. It is believed that this feature is exhibited due 
to the assumption of a constant friction factor on the rake-face. 
Uith a variable friction factor (decreasing away from the 
tool-edge), the temperature profile is likely to exhibit a single 
maximum. This aspect, however, requires further detailed 

research . 

The present study provides a theoretical methodology of 
estimating the dimensions of PSDZ and SSDZ, in a manner similar to 
that used in the analysis of photo-micrographs. The results are 
also helpful for locating the shear plane and determining the 
value of the shear angle. Further, a good estimate of the 
temperature variable near the tool-tip is obtained, with minimum 

use of empirical data. 



CHAPTER 5 


CONCLUSIONS AND SUGGESTIONS 

From the FEW predictions of the Thermal and Deformation 
phenomena in Metal Cuttinfi, the followinfi conclusions can be 

drawn : 

1. Viscoplastic Model for Metal Cutting can be successfully 
employed for an accurate prediction of the velocity and 
temperature fields in the vicinity of the tool-tip. 

2. The deformation patterns and isotherms obtained confirm the 
salient features that can expected for an Orthogonal Metal Cutting 
situation. 

3. The constant strain-rate curves plotted with cut-off values 
, 4% and 5% give a reasonably good idea of the dimensions of the 

primary and secondary shear deformation zones. By increasing this 
cut-off value, the shear plane also can be obtained which helps in 
determining the mean width of PSDZ. 

4. From the parametric study conducted, it is observed that the 
depth of cut influences the deformation and isothermal patterns 
much more than the cutting velocity or the chip-thickness ratio. 
This seems to be particularly so for the temperature variation at 
the tool-chip interface. 



5. The deformation pattern near the chip free surface seems to 
have some dependence on the shape chosen to describe the free 
surface . 

6. The choice of friction boundary condition at the chip-tool 
interface seems to have lot of influence on the temperatures on 
that boundary. A variable friction factor, decreasing away from 
the tool-edge is recommended. 

SUGGESTIONS FOR FUTURE WORK 

1. The velocity and temperature fields obtained could be further 
processed for predicting the required cutting forces which is one 
of the important input parameters. 

2. A variable friction factor could be used for the friction 
boundary condition for improving the temperature predictions at 
the chip-tool interface. 

3. Prediction of the contact length and the chip-thickness ratio 
can be attempted using the Energy Minimization principle. 

4. The effects of thermal softening and work-hardening could be 
included . 

5. The present formulation for Orthogonal cutting could be 
extended to Oblique cutting also. 
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APPENDIX A 


MATERIAL PROPERTIES 


Density (p) = 7880 

Specific heat of workpiece material = 500 J/Kg.K 

Thermal conductivity of workpiece material (k) = 50 U/mK 

Thermal conductivity of tool material (k^) = 25 U/mK 

Uniaxial Yield Stress (o ) = 300 MPa 

y 


HEAT TRANSFER CONDITION 

Overall heat transfer coefficient (h) = SO/v 

Ambient temperature ® 


CUTTING PARAMETER 


S.No. 

Cutting 

Velocity 

(m/s) 

Chip thickness 
Ratio 

Friction 

factor 

Contact 

Length 

1 . 

1.0 

0.757 

0.637 

1.85 

2. 

1.65 

0.798 

0.728 

1.85 

3. 

1.20 

0.852 

0.655 

1.47 


APPENDIX B 


DERIVATION OF EQUATION FOR VISCOSITY 


For a viacouB, incompreeaibl e fluid we can write the 
conatitutive relation in the form 

<y..=<y«5 +2#je ; 6,. = 1 i = j 

ij ij iJ ij (1) 

= 0 i J- j 

in which a is the mean stress and are the strain rates. 

Alternatively we can rewrite equation (1) as 


'ij 


2p •'ij 


C2) 


in which are the deviatoric components of stress. 

For a visco-plastic, 'associated' material we can write with 
some degree of generality 


e 


i j 


r 


dF 

Ay. . 

13 


where ^ is a constant 'pseudo-viscosity', and 


(3) 


F = = 0 


(4) 



reprcatnte th* plaetlc ylald condition. Further we use the 
following notation which ensures no plastic flow when stresses are 
below yield 

4f> = F of F > 0 and <F> =0 of F S 0 (4) 

Clearly as the constant fj ■* 0 , the equation (3) specifies the 
behaviour of an Ideally plastic mateiral. Specializing the above 
relationship to the von nises criterion for which 

F =y/(y “ J ^ 

where Y is a constant representing some yield parameter (here as 
the uniaxial test yield value) we can write 



If flow occurs, by definition (3) F iO and we can identify 
expressions (2) and (3). Using also equations (5) and (6) we can 
writ e 

( 7 ) 


®lj®i j^ 



From equation (2), we can also write the Identity 

= " ( 8 ) 

In which ela the aecond atraln Invariant and Inaertlng thla Into 
equation a.) we obtain finally an equivalent vlacoalty aa 


(i-)Y + ^^6 



e 


In which the parameter e deacrlbea an ’effective’ atraln rate. 
Relationa (9) are a particular caae of a non-Newtonian fluid where 
in general we can write 

^ = fj(e) (10) 

and where the functional form is specified in an experimental 
manner. Aa can easily be observed the viaco-plastic model is a 
generalization of the special caae ofa Bingham fluid. 
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r^r^r^oori 0000 ci o a 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


'■ ■■ i i t i i .?. ■» f .ii * 'if -ii * .-S If if ■» .f .» ? .} !i. Y :j! * * li: * !j! !» * * * Sjf * .6 * * * 

MASTER FLUID 

*.****» ^ J 1> 1 1 - : >■ t if’ r :• t • i Jf t i .i .} {f it i< :t"It i^ « it « .{t it f it 4: « if ')( fit -{i ^ a). i{: :jt $ « « * ^ iji ;(! ift $ « !(e 

3“ EXTE'.RNAL SUBRGUTIHES 

DEHENS 
DINPUT 
_ DPfjvrs 

I- ITfERAT 

C0HH0N/ARC1 /POSGPC 3 ) , WEIGP( 3) 


COMHE-N-'AFCS/LNODSC 1 05, 8> , COORD(36 6, 2 J ,NADF!1( 3S6? ,NODFH( 366 3 
C0rt?10»«/ARC3/VARB1 C 1229) ,VARB£(1EE9) 

C0HHCN/ARC35/EPSN< 105,95 
DIMENSION EQRHS{12£9) 

DIMENSION LBOUDt 1229) ,LHEDV( ISO) ,PNORt1t 150) ,BOUDV (1229) 
DIMENSION GFLUMi ISO, 150) 

OPENUJNir = £0,FILE='s3v3. Inp' ) 

OPENlUNIT = 40,FILE='var. inp‘ ) 

0PENtUNIT = 24,FILE='s3v3.out ' ) 

OPEN? UNIT *35. FILE='s3v3v. out ' 3 
0P£N(UNIT=E8,FILE='s3v3e.oUt ' ) 


*** SET DYNAMIC DIMENSION VALUES 

CALL DIHENSIMELEM, HFRON,HP01N,MTOTV3 


*** READ IN ALL PROBLEM DATA. 

CALL D INPUT! I AXSY , BOUDV , LBOUD , MELEM , 

! HPOIN,MTOTV,NBCON,NDOFH,NELEM,NEVAB,NGAUS, NITER, 

! NNODL, NNODP,NPO IN, NTOTV, RELAX, TOLER, 

! XFORC.YFORC) 

*** calculate: the shape function and derivative values 


*** SET UP EQUATIONS AND ITERATE UNTILL SOLUTIONS CONVERGE. 

CALL ITERATC IAXSY,GFLUM,BOUDV,LBOUD,LHEDV,PNORM,EQRHS, 

! HELEM,MFRON,MPOIN,MTOTV,NBCON,NDOFM, 

! NELEM.NEVAB.NGAUS, NITER, NNODL,NNODP,NPOIN,NTOTV, 

! RELAX, TOLER, XFORC,YFORC) 

URITE{£8, 100) ■ 

100 FORMATCEOX'STRAIN RATE DISTRIBUTION',//, 

! 'KELEM',4SX, 'GAUSS POINTS' ,/, 

! 18 X, ' 1 ' , 12x, '2' , 12X, '3' , lEx , '4' , 1 OX, '5' . 1 OX, '6' , 10X, '7' , 

! 12X, '8' , leX, '9' ,/) 

DO 200 I « 1 ,MELEM 

URITE(28, 1 1 1 ) I, fEPSNI I, J) , J«1 ,9) 

111 FORMAT <13, 5x ,9(E1 0 .4,2X3 3 

£00 CONTINUE 

type*, 'MAIN OVER' 

STOP 

END 

SUBROUTINE D IMENSC MELEM, MFRON, MPOIN, MTOTV ) 

C 

C*** SET DYNAMIC DIMENSION ARREY SIZES. 

C 

MELEM«105 

MFRON«t50 

MP0IN»366 



{; 4£ 1 1 1 1 

I 


I 


c 

c 

c 

01“** 


c 

£*** 

c 

c 

Cl 000 

c 

C2000 


ceoi 0 

c 

c 


EOEO 

c 

c 

c 

C«*¥ 

c 


S030 

I 

! 


|t1 1S3 

! 

! 


Sl'eftOU'^ IME 
MELEM, !1P0; 
, NSAUS^N!' 
TOLER , ;-«FC:'J 
COHHON.-'AFf 
Dins MS ION 
COfinOM 'AR< 
COJ-IMOM/ARi 
.--API 

COMMON/AP! 
C0»1M0M/AR( 
COMMON/ API 
OPEN < ON r 
OPEN( UNI' 


DINPUK 1AXSY,B0UDV,LB0UD, 
i . MTOTV , NBCON , NDOFM , NELEM , NEV AB 
:R, NMODL,NNODP|NPOIN,NTOTV, RELAX, 

YFORC) 

^'INODSC 1 02,8} , COORD (366, S) . NADFMC 366 ) , NODFH ( 366 ) 
-BOUD C MTOTV } , BOUDV( MTOTV) 

J/VARB1 ( 1££9) , VARB£( ieS9) 

1 0/CUTyEL,DCUT,BSTAR 

I £/HSTAR,TFSTAR,FRIFAC,PI ,EN, COSA,SINA,CONRAT 
I 3/USTAR,PECLET,SieMA,RHO,MBOUEL 
iA/ IBELCSO), I SURF (50), I SIDE (50) 

J0/ANSLN{366) 


iO,D£VICE = 


FILE='SMALL. INP' ) 


£4,DEVICE='DSK' , F I LE= ' SMALL . OUT ' ) 


GIVE VALUES TO THOSE VARIABLES UHICH CAN NOT BE CHANGED 

RHOs?8SO . 0 

ND0FM=4 

NEv'AB=£8 

NGAUS=3 

NNODLaA 

NNODP»S 


READ IN EACH LINE 4 ECHO IMMEDIATELY 


READ(£0, lOOOlTITLE 
FORMATS 1EA6) 

WRITE(£4,E000) TITLE 
FORMATC 1.H1 . //IX, 1EA6) 

READfEO, *)IAXSY, NELEM, NITER, NPO IN, NRPON 
pause 1 

URITE(£4, ») IAXSY,NELEM, NITER,NPOIN,NRPON 
WRl TE(E4, *1 NELEM, NPO IN, NNODP 
F0RMAT(//t3H CONTROL DATA,/13H **♦*♦♦*****,// 

8H lAXSY *,I4,4X,8H NELEM *.I4,4X,8H NITER =,I4,4X, 

8H NPOIN *,I4,4X,8H NRPON =,I4) 

READSEO, i }NI CON, NBCON, NEBCN,NNBCN 
pause E 

URIT£(24,£OaO)NICON,NBCON,NEBCN,NNBCM 
F0RMAT<8H N I CON- , 1 4 , 4X , 8H NBCON- , 1 4 , 4X , 

SH NEBCN®, I4,4X,8H NNBCN-,14) 

CHECK INITIAL DATA 

CALL DIAGN1 (NBCON, NEBCN, NELEM, NI CON, NNBCN, NPOIN, NRPON) 


READ FLOW PARAMETER 

READ ( SO, x=)RELAX, TOLER, XFORC,yFORC 
pause 3 

WRITE CS4,S030) RELAX, TOLER, XFORC, YFORC 

FORHAT(//EOH PHYSICAL PROPERTIES, /SOH *******************, 
//4X,8H RELAX - , FI 0 . 5, 4X, 8H TOLER -,F10.S, 

/8H XFORC -,F10.5,4X,8H YFORC »,F10.5) 

READ (20, * )EN,GAMMA,SIGMAY 

pause * 4 ^ 

WRITE(E4,£1 1 53 )EN, GAMMA, 81GMAY 

F0RMAT(//22H VISCOSITY PARAMETERS , /22H *********************, 
//4H EN«,F10.5,4X,7H GAMMA- , El 0 . 5 , 4X 
,8H SIGMAY“,E1 0 .5) 

READ(£0,*) CUTVEL,DCUT,RAKE,FRIFAC 



URI!£i£,, MSURF 

Fr. D SCi . » 1 ?tELS 1 . HELSa , MELS3 , MELS4 , MELSB , HELSS , HELS7 , RELSS 
pause • 1 £ •’ 

m [ IE (£4. -ts/HELSt ,i1ELS£,MELS3,M£LS4,MELS5,MELSS,MELS7,«ELS8 
f JHELS1 ,11ELS£.i1ELS3,HELS4;f1ELS5,HELS6,HELS?,MELS8 
heOU£L=H£LS1 ^HELS2+MELS3+MELS4+MELSS+MELS6+MELS7+«ELS8 
DO 9 1 1 = 1 ..HEOUEL 
n = 13 

READISO.O JJ. ISURF(in,IBEL(in,ISIDE(in 
WRnECc•4.=^) JJ, ISURFCII),IBEL(II),ISIDECII> 

WPITECS.’*) n, ISURF( ID , IBEL« U ) , ISIDEai ) 

9 COriIZNUE 

pause ' 1 3 mbouel ' 
yf>nE<£4,*) HBOUEL 
WRITE (5.0 HBOUEL 

SET UP THE VECTOR GIVING THE D.O.F. AT EACH NODE 
C 

IT£ri?=NN0DP-1 
DO 40 IELEH=*1 , NELEM 
DO 40 3N0DP*1 . ITEHP.2 
NODFfKLNODSC lELEM, INODP) )=NDOF« 

NOOFH(LNDDS( lELEH. INODP+1 J )»ND0FM-1 
40 CONTINUE 

C 

CK** INTEPPOLE mOSIDE NODE COORD WHERE NOT KNOWN 

C SKIPPING THE INTERPOLATION DUE TO AVAILABILITY OF COORDINATES 

GOTO 71 

DO 70 rELEH=1 , NELEM 
DO 60 INODP^a.NNODP.E 
IPOIN»LNODS( lELEM, INODP ) 

TEHPY«ABS( COORD { IPOIN, 1 ) )+ABS( COORD ( IPOIN, £) ) 

IFtTEMPY NE.0.0) GO TO 60 
JPOIN»LNODS{ lELEM, INODP-1 ) 

KNODP-INODP+1 
IFTKNOOPST NN0DP)KN0DP«1 
KPOIN*LNODS( lELEM. KNODP) 

DO SO IDIME*1 ,e 

COORDC IPOIN, IDIHE)=( COORD( JPOIN, IDIHEJ+COORDIKPOIN. IDIHE) J* .5 
SO CONTINUE 

60 CONTINUE 

70 CONTINUE 

71 CONTINUE 
C 

C*** SET UP VECTOR WITH THE FIRST D . 0 . F . AT EACH NODE 
C 

NADFM( 1 )»t 

DO 80 IPOIN»e,NPOIN 

NA0PH( IPOIN)»NAOFri( IPOlN-1 J+NODFM( IPOIN-1 J 
SO CONTINUE 

KT0XV»NADF!1(NP0IN)+N0DFH(NP0IN)-t 
TYPE* iNTOTV 
pauB# 'ntotv' 

C*** INITIALISE REMAING ARRAYS. 

C 

DO 90 IT0TV*1 , NTOTV 
LBOUD( 1TOTV)«0 
BOUDV(ITOTV)*0.0 
VARBl ( ITOTV)»0. 0 
VARBS( ITOTV)*0.0 
90 CONTINUE 

DO 868 ITOTV » 1 iMTOTV 
READ{40,*)VARB2( ITOTV) 

TYPE* , VARB£( ITOTV) , VARBl ( ITOTV) 

868 CONTINUE 

C 



1 00 
1 1 0 

c 

c 


c 

c** * 


S090 

I 


238 


140 
SI 00 
150 
£1 t 0 


StSO 


IfiO 

E130 

US 

8140 

170 

175 

77799 


{300 0 
;f99o 


IF« NICON . EQ . 0 UO TO 110 

DO 1 5 0 UCON=1 ,NICON 

READUO, <')1POIN, IDOFMiTEHPY 

IFUDOFM.ST 1 HDOFM=NODF«( IPOINl 

JTOTV=r4ADFM( I PO I N ) + I DOFM - 1 

VARBI ( JTOTV ) =rE!iPY 

VARBSC JTOTVUTEMPY 

CONTINUE 

CONTINUE 

pause '110' 

CHECK NODAL CONNECTION i COORDINATES 

CALL. D I ASN£{ MELEH, MP01N,NELEM .NICON, 

NNODP,NPOIN,NTOTV) 

READ IN BOUNDRY CONDITION 
pause ' bound . con ' 

WRITE{£4,S090) 

FORHAT C / •''SOH BOUNDARY CONDITIONS//EOH a******************,// 
35H NODE U-FIXED PRESSURE V-FIXED) 
pause ' reading be.' 

WRITEIS, « INBCON 
pause ' nbcon ' 

DO 175 IBCON=1 , NBCON 

READ? £0, ¥5 IPOIN, IDOFM,BVALU 

IFt IDOFH EQ . 4) BVALU = BVALU/TREF 

write (5,*) IBCON,IPOIN 

IT0TV*NADFf1( IPOIN) + IDOFM-1 

IFlNODFMt IPOIN) .EQ.4) GOTO £38 

IF( ( IDOFH+1 ) , GT .NODFH( IPOINl) ITOTV=ITOTV-l 

LBOUDI ITOTV)»1 

BOUDVI ITOTVl-BVALU 

GO TO ( 140, 150, 160, UEIIDOFM 

CONTINUE 

WRITE (£4, £100) IPO IN, BVALU 
FORHATI IS,F9 .4) 

GO TO 170 
CONTINUE 

WRITE(£4,E1 1 0 ) IPOIN, BVALU 
FORMAT! IS, 1 0X,F9 .4) 

IF(NODFHnPOIN) .NE. ( IDOFH+1 ) ) GO TO 170 

URITE(24.£1£0)IBCON, IPOIN 

FORMAT{/S7H ERROR!! BOUNDRY CONDITION, 14, 

42H DEFINES A PRESSURE AT MIDSIDE NODE NUMBER, 14) 
stop 

CONTINUE 

UR I TEC £4, £130) IPOIN, BVALU 
FORMATC I5,£0X,F9.4) 

GOTO 170 
CONTINUE 

WRITE<24,2140) IPOIN, BVALU 
FORMATC IS, 30X,F9. 4) 

JTOTV » ITOTV 

URITECS, *) JTOTV, LBOUDC ITOTV ),BOUDV( ITOTV), IPOIN 
CONTINUE 

type*, 'DINPUT OVER' 

DO 990 ITOTV » 1 ,NTOTV 
JTOTV » ITOTV 

UR I TEC £4, 3000 ) JTOTV, LBOUDC ITOTV) , BOUDVC ITOTV ) 

FORMATC lx.' IT0TV=' , I3,£x, 'LBOUDC ITOTV)*' , I3,£x, 

'BOUDVC ITOTV)*' ,FtO. 4) 

CONTINUE 

DO 300 1*1 ,MPOIN 



300 


350 


400 




c 

c*** 

10 

c 

c 


C 0 ri i I N U E 

n = ISO 

J 1 = £ i 4 

LI ® 14-^ 

DO 3E0 1=1,6 

1 1 = n 1 

Ji = jt + 1 

11 = LI + 1 

ANGL1H I ! i * £0.0 

AMCLf'H JI > »■ £0.0 
AHGL»HL 11 » 90.0 
CONTINUE 

ANGINISS! = 90.0 
ANGLN1S3i = 90.0 
ANGLN(e4> = 90.0 
ANGLM 144) = 90.0 
ANGLNf £53) =85.0 
ANSLNC£54} = 84.0 
ANGLN1255) = 77.0 
AN&LN<£56) =r 75.0 
ANeLN(£57) = 68.5 
ANeiNlSBS) = 67.5 
ANGLN<£59) = 66.25 
ANGLN1 £60) =65.0 
AN&LN1S61 ) - 63.0 
ANGLN126S) = 60.5 
ANGLN<£63) * 54.0 
ANGLN(264) = 49.0 
ANGLN{271 ) = 20.0 
ANGLN{£7£) = 20.0 
ANGLW(337) » 45.5 
ANGLN1338) » 42.0 
ANGLNI 339 ) * 41.0 
ANGLN(340) ■ 37.0 
ANGL.N( 341 ) » 30,5 
ANGLN(342) » £0.0 
ANGLN(343) « £0.0 
ANCLN(344) * £0,0 
ANGLN(363) » 90.0 
ANGLN{364) « 90.0 
ANGLN(365) * 90.0 
ANGLN(366) ■90.0 
DO 400 KKL « 1 ,HPOIN 
LL » KKL 

Typ*»,LL, ANGLN(KKL) 

CONTINUE 

RETURN 

END 

SUBROUTINE D I AGNK NBCON , NEBCN, NELEH, NICON , NNBCN, NPOIN, NRPON ) 
DIMENSION NECHO( 150) ,NER0R(3) 


SCRUTINY OF CONTROL 
DO 10 IEROR»1,8 
NERORI IEROR)«0 


DATA 


SCRUTINISE CONTROL DATA AND PRINT ERROR MESSAGE, 


IF(NPOIN.LE.O) NEROR(n = 1 
IF(NRPON.LE.O,OR.NRPON.GT.NPOIN)NEROR<£)=1 
IF( NELEM . LE . 0 ) NER0Rt3)®1 
IF1NP0IN.GT.NELEM*8) NER0R(4)=1 
IF(NBC0N.GT.NP0IN*4) NER0R(5)=1 
IF{NIC0N.GT.NP0IN*4) NER0R<6)*1 
IFCNEBCN.GT. NELEM) NER0RC7)»1 


c 

c* * * 
c 


c 

c 

2000 

ao 

c 

c 

aoi 0 

30 

1 000 

aoao 

c 


Q lie :|i 41 * 


C 

c*** 

c 

c# 4> 4i 

c 


1 0 
c 

c** * 
c 
c 
c 

c*«« 

c*** 

c 


ao 

30 

40 

c 

c*** 

c*** 

c**« 

c 


CHECi< m ERROR-IF UNITY PRINT. 

JEROR=0 

DO as IE'P0R=1 ,8 
!F«NERORriEROR) ,EQ. 0)GO TO 20 
JER0R=1 

WRITE 0U1 ERROR NUMBER. 

UR!TE(a4, iiOIEROR 

FORNiST t /, 1 0X,33H CONTROL PARAMETER ERROR****ERROR , IS ) 
CONTINUE 

IF( JEROR .EQ. ORETURN 

LIST DATA REMAINING AFTER CONTROL PARAMETERS. 

WRIT£{a4, «“) 

FORMAT (/J OX, 37HDATA FOLLOWING ERROR IN CONTROL CARDS/) 
CONTINUE 

READ! £0, 1 000 )NECHO 
FORMAT! 1E0A1 ) 

WRITE(£4, *)NECHO 
FORMAT! 1 OX, 80A1 ) 

GO TO 30 

type* , 'DIAGl OVER' 

RETURN 

END 

SUBROUTINE D I AGNE ! MELEM , HPO IN , NELEH , N I CON 
,NNODP,NPOIN,NTOTV) 

DIMENSION NECHOt 150 ) ,NEROR( 13) 

COMMON/ARCa/LNODS( 1 05,8) , COORD ( 366 , E ), NADFM C 366 ) ,NODFM{ 366) 

SCRUTINISE NODAL POINT AND ELEMENT DATA. 

INITIALISE ERROR ARRAY. 

DO 1 0 IEROR=l , 13 
NERORt lEROR )*0 
CONTINUE 

CHECK PHYSICAL PROPERTIES. 

IF { DENS . LE . 0 . 0 . OR . V ISCY . LE . 0 . 0 ) NEROR ( 9 ) =1 

CHECK NODAL COORDINATES. 

CHECK FOR TWO IDENTICAL COORDINATES. 

DO 40 IPOIN»e,NPOIN 

JPOIN*IPOIN-1 

DO 30 KPOIN-1 , JPOIN 

IF { COORD! IPOIN, t ) . NE . COORD ( KPOIN, I ) . OR . COORD ( IPO IN, a ) . NE . 

COORD(KPOIN,a) )GO TO £0 

NEROR(10)«1 

CONTINUE 

CONTINUE 

CONTINUE 

CHECK ELEMENT NODAL NUMBERS 

A — REPITITION OF ANODE NUMBERIN ONE ELEMENT 
B— -NOD NUMBER OUTSIDE PERMISIBLE BOUNDS. 

DO SO IPOIN»1 ,NPOIN 
DO 50 IELEM=1 ,NELEM 



DO £0 INC'DE'-I ,, NNODP 

IF5tNODSf lELEn, INC1DE5 .NE. IPOIN)GO TO 50 
LC0NT=LC0NT+1 
IFCLCOMT , GT . 1 JNERORC 1 1 ) = 1 
so CONTINUE 

C GOTO GO 

DO £0 IELEM«1 , NELEM 

DO SO 1NCDE=] ,NNODP 

IF) LNODSC IElEM, INODE ) . LE . 0 . OR . LNODSC lELEM, INODE ) . GT . NPOINJ 
! NEROR<1£.* = 1 
60 CONTINUE 

C 

C*** CHECK ON NUMBER OF INITIAL CONDITIONS 

C 

IF(N1C0N.5T. NTOTVINERORC 13)=1 
C 

C*** WRITE OUT ERROR NUMBER 

C 

JEROR=0 

DO 70 lEROR'lO, 13 

IFINERORC lEROR) .EQ. 0)GO TO 70 

JEROR«1 

WRITE{£4, ^IIEROR 

EOOO FORMAT1 /I OX, 1 OHDATA ERROR, 15) 

70 CONTINUE 

C 

DECIDE WHEATHER TO CONTINUE OR ECHO DATA. 

C 

IFIJEROR EQ.O)RETURN 
C 

C*** LIST REMAINING DATA 

C 

URITECe4, f) 

£010 FORMAT! /I OX, 14HREMAINING DATA/) 

80 CONTINUE 

READ(£0, 1 000)NECHO 
1 000 FORMAT! 150A1 ) 

WRITE(S4,*)ECH0 
£020 FORMATC/1 0X,80A1 ) 

C GO TO BO 

type*, 'DIAGN£ OVER' 

RETURN 

END 

SUBROUTINE I TERAT ( I AXSY , GFLUH , BOUDV , LBOUD , LHEDV , PNORH , EQRHS 
! ,MEL£M,MFRON,!1POIN,MTOTV,NBCON,NDOFM, 

! NELEM, NEVAB, NGAUS,NI TER, NN0DL,NN0DP,NP0IN,NT0TV, 

! RELAX, TOLER, XFORC,YFORC) 

c__ 

cl EXTERNAL SUBROUTINE 

C_ PRESCR 

C_ FRONTS 

C_ WRITER 

C_ TOLREL 

C_ 

C**************************** 

C0MM0N/ARC1/P0SGPC3) , WEIGP(3) 

COMMON/ ARCS/LNODS! 105,8), COORD! 366, S) , NADFH(366) ,NODFM{366) 
COMMON/ARC3/VARB1 ( 1 £29 ) , VARB2 ( 1 229 ) 

C0MH0N/ARC35/EPSN{ 105,9) 

DIMENSION EORHSCMTOTV) 

DIMENSION LBOUD {MTOTV),LHEDV(MFRON),PNORMtMFRON),BOUDy{MTOTV) 
DIMENSION GFLUH!HFRON,HFRON) 

C 


c 

I ITER^O 

DO 1 1 = 1 . riELEM 
DO 1 J »■ 1 , 9 
EPSNC I , J3 = 0.0 

1 CONTINUE 

10 CONTINUE 

1 1 TE!?=I I TE!? + 1 
TYPE*, ' I2TEP=' , IITER 
WRITEOS.*)' IITER=S IITER 

DO 20 IT0TV=1 ,NTOTV 
EQRHSC ITOTV)=0 . 0 
20 CONTINUE 

C 

C*** CALL FRONTS TO SET UP AND SOLVE GOVERNING EQUATIONS 

C 

CS7 F0R!1AT(4X, ' YES 4'J 

DO 50 1=1 , NFRON 
DO SO J=1 , MFRON 
GFLUnC I , J J = 0 . 0 

so CONTINUE 

KOUNT = 1 
TYPE»',KOUNT 
NESPBC * 0 

CALL FRONT S? lAXSY, 1 1 TER , HELEN , I1FRON , NPOIN , fITOTV , NBCON , NELEM , NEVA 
! B, NNODL, NNODP , NPO I N , NTOTV , XFORC , YFORC , NGAUS , NESPBC , GFLUM , 

! BOUDV, LBOUD. LHEDV,PNORH,EQRHSJ 
KOUNT = KOUNT + 1 
TYPE*, KOUNT 

C*** CALL WRITER TO OUTPUT ITERATION RESULTS, 

C 

c CALL WRITER? IITER, HPOIN, NTOTV, NDOFM,NPOIN, HELEN! 

C 

KOUNT « KOUNT + t 
TYPE*, KOUNT 

C*** CALL TOLREL. TO CHECK CONVERGENCE AND RELAX VALUE IF NOT. 

C 

CALL TOLREL? 1 1 TER , NTOTV , NCONV , NTOTV , RELAX , TOLER , VARB 1 ,VARB2 
! ,NADFN,NODFM,MPOIN) 

C 

KOUNT = KOUNT + 1 
TYPE*, KOUNT 

C*** RETURN TO MASTER NO. IF ITERATIONS EXCEED MAXIHUM. 

C 

IFCNCONV EQ .1)60 TO 12 
IF? IITER ,LT.NITER)GO TO 10 
WRITECS, 20005 
c WRITE(35,2000 ) 

2000 FORMAT?//3£H SOLUTION HAS FAILED TO CONVERGE) 

12 CALL WRITER? IITER, HPOIN,MTOTV,NDOFM,NPOIN,MELEM) 

type*, 'ITERAT OVER' 

RETURN 

END 

SUBROUTINE WR I TER ? 1 1 TER , MPOIN , MTOTV , NDOFM , NPOIN, HELEN) 
C0MM0N/ARC2/LN0D8C1 0S,8) , COORD? 366, 2 ) , NADFM? 3Sfi ) , NODFHC 3fi6 ) 
C0NH0N/ARC3/VARB1 (1 £29 ) , VARBE? 1 £29 ) 

C OPENCUNIT»3S,FILE«'8HALL1 .OUT' ) 

«RITE(5, IIDIITER 
111 FORMAT?1X, ' IITER*' , IS) 

WRITE(5,£000)nTER 

WRITE?3S,£000)UTER 

£000 F0RMAT{//£9H RESULTS FOR ITERATION NUMBER, I£, 

• //14X,4?3HNEU,7X),7X,4C3H0LD,7X)/, . 

! 7H NODE , E< 4X , 42HU-VEL0CITY PRESSURE V-VELOCITY TEMPERATURE)) 



IODFM=NO&FH^ IPOIN J 
IADFH-NADFi1< IPOIN) 

IFUODFM. EQ , NDOFH)SO TO 10 

WRITEC E, £01 0 HPOIM, ( VARBt ( I ADFM + JODFH- ! ) , JODFM=1 , lODFH) , 

I {VARB£( lADFM+JODFM-l ) , JODFM=1 , lODFM) 

WRITB13E. £010)1POIN, ( VARB1 ( I ADFM+JODFH-1 ) , JODFH=1 , lODFH) , 

! CVARESC IADF!14 JODFH-1 > , JODFH»1 , lODFM) 

GO TO £0 
10 CONTINUE 

WRITE C 5, e0£0 3 IPOIN, ( VARBl ( I ADFH+ JODFH- 1 ) , JODFH=l , lODFH) , 

! {VARBEt I ADFH+ JODFH-1 ) , JODFH=1 , lODFM) 

WRITE 135, £CI£0)IPOIN, 1 VARBl C I ADFH+ JODFH-1 ), JODFH® 1 ,IODFH), 

« {VAR8£( IADFM+JODFH-1 JODFM=1 , lODFH) 

so CONTINUE 

SOI 0 FORMAT { I S , 4X , S C E 1 0 . 3 , 1 OX , £E 1 0 . 3 , 6X J ) 

SOHO FORMATUfe,4X,S{4E10,3,6X3 ) 

RETURN 

END 

C FINISH* 

SUBROUTINE TOL REL ( 1 1 TER , HTOTV , NCONV , NTOTV , RELAX , TOLER , 
i VARBl , VARB£,NADFH,NODFH,HPOIN) 

DIMENSION VARBl ( MTOTV ) , VARBS ( MTOTV } , NADFM { MPO IN 1 , NODFM ( MPO I N ) 
C 

C*** CHECK TO SEE IF SOLUTION HAVE CONVERGED. 

C 

C*** CHECK FOR CONVERGENCE TO REQUIRED TOLERANCE. 

C 

NDOFM ■ 4 

DO SOS IPOIN«1 , MPOIN 
IODFM»NODFM( IPOIN) 

IADFH»NADFM( IPOINJ 

IF( lODFM . EQ . NDOFM )GO TO 108 

WRITE! 5, 301 0 ) IPOIN, C VARBl ( I ADFM+ JODFM-1 1 , JODFM»i , lODFM) , 

! IVARBSl IADFM+JODFM~1 ), JODFH®! , lODFH) 

GO TO SOS 
108 CONTINUE 

WRITE! 5, 30S0 ) IPOIN, ! VARBl ! I ADFH+ JODFH-1 ), JODFH® 1 , lODFHl, 

! { VARBS! I ADFn+ JODFH-1 3 , JODFH® 1 , lODFHJ 

208 CONTINUE 

3010 FORMAT! 1 6 , 4X , 2 ! E 1 0 . 3 , 1 OX , SE 1 0 . 3 , 6X ) ) 

30E0 FORMAT! 16, 4X, S! 4E1 0 . 3, 6X) J 
CANLA®0 . 0 
NCONV®1 

UMAX » VARBl < 1 } 

PMAX « VARBl CS) 

VMAX ® VARBl (3) 

TMAX » VARBl !4) 

DO 10 IPOIN®! , MPOIN 
lODFH » NODFMdPOIN) 
lADFH « NADFM! IPOIN? 

IF! IO0FH . EQ . 3 J GOTO 16 

IF! ABSIVARBI ! lADFH) } .GT.ABS(UMAX) ) UMAX = VARBHIADFM) 

IF! AB5CVARB1 ! IADFH+1 } ) -GT. ABS(PHAX> J PHAX * VARBl ( 1ADFH+ 1 J 
IF! ABS! VARBl ( lADFH+EJ J . GT.ABS(VHAX3 J VHAX = VARBl ( I ADFM+H J 
IF!AeS(VARBl!IADFH+3J).GT.ABS{THAX)J TMAX « VARBU IADFH+3) 
GOTO 100 

16 IFCABSCVARBI ! lADFMJ J .GT, ABStUHAX) J UMAX “ VARBltlADFHJ 

IF! AB5(VARB1 ! lADFH+l ) J .GT. ABS(VHAX3 ) VHAX * VARB 1 ( IADFM+ 1 J 
IF! AB5! VARBl! lADFH+E ) ? . GT . ABS! THAX3 3 THAX ® VARBl ( IADFM+S3 
100 CONTINUE 

t yp©* , UMAX , PHAX , VHAX , THAX 
10 CONTINUE 

t yp©* , UMAX , PHAX , VMAX , THAX 


19 

34 


fiO 


70 


80 


18 


aooo 

! 


£01 0 


! 


20 

C 

C*»* 


FCUf - e V*P!1AX 

VCUT = 0 1^V«AX 

TCU7 * 0 1*THAX 

tYpe^ .UCUT.PCUT, VCUT, TOUT 

DO IS IPOIN = 1 ,MPOIN 

lODFfI - NODFM(IPOIN) 

lADFSI = NAOFHC IPOIN) 

ITOru = EADFtI 
IF( lOOFH EQ . 3) GOTO 19 
ITOTP » IADFM+1 
ITGTV » ITOrU + £ 

ITOTr * ITOTU + 3 
GOTO 34 

ITOTV * ITOTU +1 
ITOTT » ITOTU +£ 

IF(AeSi VAR81 C ITOTU) ) .LT.ABSCUCUT) ) GOTO SO 

IF( ABS( VARBI C ITOTU) ) .LT . 0 . 01 )GO TO 60 

CANGE*ABS( < VARB1 ( I TOTU ) -VARB2 { I TOTU ) J/VARBI { ITOTU) ) 

IFtCANGE GT . TOLER)NCONV«0 

IFICANLA GT.CANGE)GO TO 60 

CANLA»CANGE 

LTOLA=ITOTU 

CONTINUE 

IFtlODFM EQ.3) GOTO 70 

IF( ABSi, VARB1 ( ITOTP )) .LT .ABSCPCUT) ) GOTO 70 

IF(ABS(VARB1 ( ITOTP) ) .LT.O. 01 )GO TO 70 

CANGE*A8£( { VARB1 { I TOTP ) ~VARB2 ( I TOTP ) ) /VARB 1 (ITOTP) ) 

IFtCANGE GT . TOLER )NCONV»0 

IF(CANLA GT,CANGE)CO TO 70 

CANLAaCANGE 

LTOLAalTOTP 

CONTINUE 

IF( ABS< VARB1 ( ITOTV) ) .LT .ABSCVCUT) ) GOTO 80 

IF( A8S( VARB1 t ITOTV) ). LT . 0 . 01 )GO TO 80 

CANGE -ABSI 1 VARB1 C I TOTV ) -VARB2 ( ITOTV) )/VARB1 (ITOTV) ) 

IF(CANGE GT TOLER ) NCONV«0 

IF(CANLA GT.CANGE)GO TO 80 

CANLA-CANGE 

LTOLA»ITOTV 

COWT INUE 

IF(ABS(VARB1 ( ITOTT) ) .LT.ABS(TCUT) ) GOTO 18 

IF( ABSI VARB1 ( ITOTT) ) .LT . 0 . 01 )GO TO 18 

CANGE»AeS( { VARB1 ( I TOTT ) --VARB2 ( I TOTT ) ) /VARB1 (ITOTT) ) 

IF(CANGE GT, TOLER )NCONV*«0 

IF{ CANLA . GT .CANGE)60 TO 18 

CANLA-CANCE 

LTOLA-ITOTT 

CONTINUE 

typ#*, '5000' 

WRITE(S,2000)LTOLA,CANLA 
WRITE(35,£0C»0)LTOLA. CANLA 

F0RMAT(/37X, ' LARGEST CHANGE OCCURE ATDOF NO ,15,/, 

10X,' CHANCE »',F10.4) 

IF(NCONV.EQ.O)GO TO 20 
t yp#* , ' 6000 ' 

URITE(S,201 0) IITER 

HRITE(35,S010) IITER rni coAMrF/ 

F0RMAT<//4SH SOLUTION HAS CONVERGED TO REQUIRED TOLERANCE/, 

3x, 'IN', 3x, IS, 3x. 'ITERATIONS') 

typ#*, ‘T'OLREL OVER' 

RETURN 

CONTINUE 

IFCIITER.EQ.I )GO TO 40 

RELAX VARIABLES FOR NEXT ITERATIONS 



, t vpB-f , ' 7000 ' 

DO ESf? IP0IN=1 ,MPOIN 
I0DFf1=*ls}0DFM( IPOIN) 
lAOFMaNADFMt IPOIN) 

DO ai JODFM=1 , lODFH 
TEMPYI-VARBI ( I ADFM+ JODFM- 1 ) 

TE!1PY£*VARBe( I ADFM + JODFM- 1 ) 

VAR&£( I ADFHt-JODFH ~1 ) =TEMPYe+RELAX* < TEHPY 1 -TEMPYS ) 

£1 CONTINUE 

SEE CONTINUE 

GOTO ££5 
40 CONTINUE 

DO 50 1T0TV«1 ,NTOTV 
VARBSt !TOTV)»VARB1 ( ITOTV) 
so CONTINUE 

EES CONTINUE 

RETURN 
END 

SUBROUTINE FRONTS! I AXSY, 1 ITER , HELEN, MFRON, MPOIN, MTOTV, NBCON, NE 
! LEH , NE VAB , NNODL , NNODP , NPOIN , NTOTV , XFORC , YFORC , NGAUS , NESPBC, 6FLU 
! H, 80U0V, LBOUD,LHEDV,PNORM,EQRHS) 

c_ 

C EXTERNAL SUBROUTINES 

cl MATRIX 

C_ 

DIMENSION LBOUDIHTOTV) , LHED V ( MFRON ) , PNORM ( MFRON ) , BOUDV { MTOTV ) 
DIMENSION EQRH5!HT0TV> 

COMMON/ ARC2/LN0D5! 1 05,63, COORDt 366, 2 3 , NADFM ( 366 3 , NODFH { 366 3 
COMMON/ ARC3/V ARBI ( 1 229 3 , VARB2 C 1 ESS J 

DIMENSION LOCEL<253 ,NDCST<283 , FLUMX C 28 , 26 3 , CFLUM ( MFRON , MFRON 3 
OPEN (UN I T«£S, FORM’S 'UNFORMATTED ' , F I LE»» ' TTRASH . DAT' 3 
IF{ IITER . ST. 1 360 TO 40 
C 

C*** ON FIRST ITERATION ONLY FIND LAST APPEARANCE OF EACH NOD. 

C 

DO 30 IPOIN-1 , NPOIN 
LASTE»0 

DO 20 IELEM»1 , NELEM 
DO i 0 INODP-1 , NNODP 

IF(LNODS( lELEM, INODP) .NE. IPOIN3GO TO 10 
LASTE»IELEH 
LASTN»INODP 
60 TO 20 
10 CONTINUE 

20 CONTINUE 

LH0D8C LA5TE, LASTN3»-1P0IN 
30 CONTINUE 

40 CONTINUE 

C 

C*** INITIALISE HEADING AND GRAND FLUID MATRIX. 

C 

C REWIND 25 

C REWIND 10 

NCRIT-MFRON-NEVAB 

NFRON-0 

DO SO IFRON 1 , MFRON 
DO 50 JFRON » I , MFRON 
GFLUMC IFRON, JFRON3 « 0.0 
50 CONTINUE 

KELEM»0 
C 



^ 1.M ruwnxHto c.i-cncNim_ nwirtivco 

C 

60 CONTINUE 

KEl,EM»KELEH+l 

CALL NATR IX flAXSY, HELEN, EQRHS, FLUNK, HPO IN 
! .MTOTV.NEVAB, NNODL , NNODP , XFORC, YFORC 
! , KELEH, NELEN, NGAUS, I ITER} 

KEVAB^O 

C 

C*** CRFAI GLOBAL DOF ARRAY FOR EACH LOCAL ELEMENT DOF. 

C 

DO 70 INODP*1 , NNODP 
KPOIN^LNODS? KELEH, INODP J 
IADFM»NADFM( IABS(KPOINJ ) 

LODFM*NODFMf lABSCKPOIN) ) 

DO TO I ODFM»1 , LODFM 
KEVAB-FEVAB+1 

LOCELf KEVAB)=IADFH+IODFH~l 

IF<KFOIN LT OJLOCEL(KEVAB}»-LOCEL(KEVAB) 

70 CONTINUE 

C 

C*** FIT EACH DOF INTO THE FRONT WIDTH EXTENDING IF NECESSARY. 

C 

DO 1 £0 IEVAB*1 , NEVAB 
KTOTV=LOCeLC lEVAB) 

IFCNFRON.EO 0)GO TO 90 
DO 80 IFRON»1 , NFRON 
KFR0N»1FR0N 

IF( IABS( KTOTV} . EQ. IABS(LHEDV(KFRON} } }G0 TO 110 
SO CONTINUE 

90 CONTINUE 

NFRON-NFRON+1 

IF 5 NFRON . LE . MFRONIGO TO 100 
WRITE(5, EOOOJ 

£000 FORMAT(//3X, 'PROGRAM HALTED FRONTUIDTH IS TOO SMALL') 

STOP 

100 CONTINUE 

NDESK lEVABl-NFRON 
LHEDV(NFRON)*KTOTV 
CO TO 120 
110 CONTINUE 

NDEST< IEVAB)»KFRON 
LHEDV? KFRON )»KTOTV 
120 CONTINUE 

C 

C*** ASSEMBLE NEW ELEMENT INTO GRAND FLUID MATRIX. 

C 

DO 130 IEVAB-1 , NEVAB 
IFRON*NDEST(IEVAB) 

DO 130 JEVAB*1 , NEVAB 
JFRON*NDEST( JEVAB3 

CFLUMt JFRON, I FRON 3 *CFLUH ( -IFRON , IFRON } +FLUHX( JEVAB , lEVAB) 

130 CONTINUE 

IF(NFRON.LT.NCRIT.AND.KELEM.LT.NELEM)GO TO 60 
140 CONTINUE 

NF8UM-0 
PIVOT»0 . 0 

C*** CHECK LAST APPEARANCE OF EACH DOF PROCESS BOUMDRY CONDITIONS. 

C 

DO 170 IFRON-1 ,NFRON 

IF(LHEDVUFRON) .GE.03GO TO 170 

NFSUM-1 ^ 

IF{LB0UDUAB8(LHEDV(IFR0N))).NE. nSO TO 160 
KT0TV»IABSCLHEDV<IFR0N) } 

LBOUDCKTOTVJ — l 


150 

160 

C 

C*** 

c 


170 


C 


E010 


180 

C 

c*** 

c 


190 


C 

c 


EOO 

E10 


SSO 

S30 


E-»0 

E50 


E60 

E70 


E80 


E90 


DC 15 0 L=-R0N==1 .NFPON 
GFLUH< IFRON .LFRON5=0 . 0 
CONTINUE 

GFLUHI IFRON, IFRON)=1 . 0 
CONTINUE 

SEARCH FOR LARGEST PIVOTAL VALUE 

PIVOG=GFLUM( IFRON, IFRON) 

IF< AES(PIVOG) , LT. ABS(PIVOT) )GO TO 170 

PIVOT=PIVOe 

LP1VT=IFR0N 

CONTINUE 

IFtNFSUH EC! 0)GO TO 60 
KT0TV*1ABS(LHEDV(LPIVT) ) 

REDUCED THE VALUE IE-7 TO IE-10 
IF! ABS( PI VOT » , GT 1.0E~1£)60 TO 180 
TYPE+ , PIVOT 

WRITEIS, E0101KTOTV, PIVOT 

FORMAT t //3SH PROGRAM HALTED I LL-COND I T I ON 1 N6 , // 1 7H D.O. FREEDOM 
.I4./13H PIVOT VALUE ,E9.E) 

STOP 

continue 

NORMALISE PIVOTAL EQUQTION 
DO 1 90 IFRON*l .NFRON 

PNORMt IFRON5»GFLUM(LPIVT, IFRON) /PIVOT 
CONTINUE 

RH5I0»EQRHS( KTOTV)/PIVOT 
EQRMS{KT0TV)*RH5ID 

ELEMINATION OF PIVOTAL EQUATION REDUCING FRONT WIDTH. 

1F< LPIVT , EO UGO TO 250 
DO E40 IFRON- 1 , LPIVT-1 
FACOR«GFLUM{ IFRON, LPIVT) 

IFiFACOR EQ.OJGO TO 210 
DO 200 JFRON- I , LPIVT-1 

GFLUM( IFRON, JFRON)»GFLUM( IFRON, JFRON ) -FA COR*PNORM{JFRON ) 

CONTINUE 

CONTINUE 

IFILPIVT EQ.NFRONJCO TO 230 
DO 220 JFRON-LPIVT+1 , NFRON 

GFLUM( IFRON, JFRON-1 )*GFLUH( IFRON , JFRON ) -FACOR*PNORM ( JFRON ) 

CONTINUE 

CONTINUE 

IT0TV«IAB5<LHEDV( IFRON) ) 

EQRHSt 1T0TV)»EQRHS{ I TOTV ) -FAC0R*RH8 1 D 

CONTINUE 

CONTINUE 

IF(LPIVT.EQ.NFRON)GO TO 300 
DO 290 lFRON-LPIVT+1 , NFRON 
FAC0R»6FLUM( IFRON, LPIVT) 

IFILPIVT . EQ. 1 )GO TO 270 
DO 260 JFRON*! ,LPIVT-1 

GFLUMC IFRON-t , JFRON ) «CFLUM{ IFRON, JFRON)-FACOR*PNORMIJFRON) • 

CONTINUE 

CONTINUE 

DO 280 JFRON-LPIVT+1 , NFRON 

GFLUMC IFRON-1 , JFRON-1 )«SFLUM( IFRON, JFRON)-FACOR*PNORH(JFRON) 
CONTINUE 

IT0TV»1ABS(LHEDVC IFRON) ) 

EQRMSC ITOTVJ-EQRHSC 1T0TV)-FAC0R*RHSID 

CONTINUE 


c 

c*»* WRITE 0U1 r40NFIXED PIVOTAL EQUATION ON TAPE. 

C 

IF{LeO!JO<KTOTV) .NE. 0)80 TO 3t0 

URITE(E5)NFR0N,LPIVT, ILHEDVC IFRON) ,PNORM( IFRON) , IFRON=l ,NFRON) 

310 CONTINUE 

DO 320 IFROr4 = 1 ,NFRON 
GFLUN? IFRON, NFRON)«0 . 0 
GFLUIK NFRON, IFRON) = 0. 0 
3£0 CONTINUE 

IF1LPIVT EQ.NFRONISO TO 340 
DO 330 IFRON=LPIVT,NFRON~1 
LHEDVC IFRON)=LHEDVC IFRON+t ) 

330 CONTINUE 

340 CONTINUE 

NFRONwNFRON-1 

C 

C*** ASSEMBLE ELIMINATE OR BACK SUBSTITUTE. 

C 

IF(NFRON,ST.NCRIT)GO TO 140 
IFCKELEM LT.NELEfDGO TO 60 
IFtNFRON.GT 0)GO TO 140 
C 

CC»** BACK SUBSTITUTION 
C 

DO 350 ITOTV=1 ,NTOTV 
VARB1 { ITOTV J*BOUDVI ITOTV) 

LBOUDC ITOTV)=-LBOUD( ITOTV) 

350 CONTINUE 

DO 370 ITOTV»1 ,NTOTV~NBCON 
BACKSPACE 25 

READ(25)NFRQN, LPIVT, (LHEOVI IFRON ) , PNORM { I FRON ) , IFRON=1 ,NFRON) 
KTOTV»IABS<LHEDV(LPIVT) ) 

TEf1PR*0 . 0 
PNORM(LP1VT)»0.0 
DO 360 IFR0N»1 , NFRON 

TEMPR-TEMPR-PNORMC IFRON)*VARB1 C IABS(LHEDV( IFRON) ) ) 

360 CONTINUE 

VARB1 (KTOTV)«EQRHSCKTOTV)+TEMPR 
BACKSPACE as 
370 CONTINUE 

type*, 'FFRONTS OVER' 

RETURN 

END 

C)|e9|i«:i|e¥)|s«i)|<<('4:!|:%*t‘%>t<*>(i|:%¥¥4'>li9‘’|i’i<*<|:V’l‘***************************************** 

Cle,|e:|c9i4t,fi«iliyi,|i,|(if;,|tJt:4C)fiV¥V’(’l<*T’l")‘****************************************** 

SUBROUTINE MATRIX! lAXSY, MELEM, EQRHS , FLUMX, 

! MPOIN, MTOTV, NEVAB, NNODL, NNODP , XFORC 
! , YF0RC,KELEM,NELEM,N6AUS, IITER) 

C *#*INSERTED THE DUMMY VARIABLE "NSIDE" ABOVE******** 

DIMENSION CARTL<£,4) , CARTPI £ . 8 ) , ERHSU! 8 ) , ERHSV( 8 ) , FLUMX! £8 , £8 ) 
DIMENSION ERHSTtS) , VISCIS) 

DIMENSION SHAPLI4) ,SHAPP(8) ,AREAW(9) 

COMMON/ARCe/LNODS 1105,8), COORD! 366 , £ ) , NADFM! 366 ) , NODFM! 366 ) 
DIMENSION EQRHSIHTOTV) 

C0MM0N/ARC3/VARB1 ( 1££9) ,VARB£( 1££9) 

C0MM0N/AREA3/CARPG ! a, 7e ) , SHAPG! 7£ > , CARLS I £, 36 ) , SHALG! 36 ) 

C C0MM0N/ARC4/AREAW<9) 

COMMON/ARC1 0/CUTVEL , DCUT, BSTAR 

COMMON/ARC1 a/HSTAR , TFSTAR, FRIFAC, PI , EN, COSA , SINA, CONRAT 

COMMON/ ARC 13/USTAR,PECLET, SIGMA, RHO,MBOUEL 

C0MM0N/ARC14/ IBELI50), ISURFI50), ISIDE!50> 
COMMON/ARC35/EP8N(105,9) 

C 0PEN!UNIT«a8,DEVICE»'DSK',FILE»'TRASH.0UT' ) 

type*, (LNODS!KELEM,IJK) , IJK«1 .8) 


onn 4^ w ooortf -* ooo 


% f VI* * « X" i' 4 . W »! "f,! Kl fc 1 A ^ 

CALL DRIVES^ CCS0RD,LN0DS,MELEM,MP01N,NELEH,3,NN0DL,NN0DP, 
! KELEH,AREAy) 

** INITIALISE ARRAYS 

DO to IN0DP=1 ,NNODP 
ERHSUC INQDP ) =0 . 0 
ERHSV( INODP )»0 . 0 
ERHSK INODP )®0 0 
CONTINUE 

DO £0 IEVA8=1 ,NEVAB 
DO SO JEVAB=t , NEVAB 
FLUNXC lEVAB, JEVAB)*=0 . 0 
0 CONTINUE 

« « 4! 4 c « * * t « It t » •I' 4 ! r »1 If *1 * .*■»'***«** * ******** * 

*** LOOP TO CARRY OUT GAUSS INTEGRATIONS 

DO 100 IGAUS«1,S 

KI6S!=IGAUS 

DAREA=AREAW(KIGS) 

DO 30 INODP=»1 , NNODP 

SHAPP( INODP)«=SHAPG(NNODP*< IGAUS-1 J+INODP) 

DO 30 IDIME=1 I S 

CAPTPC IDIME. INODP)= C ARPG ( I D I ME , NNODP * ( I GAUS- 1 5 + INODP) 

0 CONTINUE 

DO 40 INODL=1 ,NNODL 

SHAPL( INODL JsSHALGCNNODL*! IGAUS-1 ) + INODL) 

DO 40 IDIPIE»1,2 

CARTLt I DIME, I NODL J «CARLG C I D I ME . NNODL* C I GAUS- 1 ) + INODL) 

0 CONTINUE 

IELEM*KELEM 
UVELY=0 , 0 
RADUS»0 . 0 
VVELY*0 . 0 

*** EVALUATE RADIUS AND PREVIOUS VELOCITIES AT GAUSS POINTS 

DUDX » 0.0 

DUDY « 0.0 

DVDX » 0.0 

OVDY "0.0 

DO SO INODP-1 , NNODP 

KP0IN*IAB5( LNODSIKELEHi INODP J ) 

ITOTU«NADFM{KPOIN) 

ITOTV«ITOTU+NODFM(KPOIN )-S 
SHAPE«SHAPP( INODP) 

UVELY=UVELY+VARB£C1T0TU)*SHAPE 
RADUS"RADUS+COORD{KPOIN,E)*SHAPE 
VVELY«VVELY+VARBE( IT0TV)*5HAPE 
CARXI»CARTPt 1 , INODP J 
CARYI"CARTPCS, INODP) 

UVEL » VARB2IIT0TU) 

VVEL ■ VARBEUTOTVJ 
DUDX “ DUDX + CARXI*UVEL 
DUDY « DUDY + CARYI#UVEL 
0VDX » DVDX + CARXI*VVEL 
DVDY » DVDY + CARYI*VVEL 
50 CONTINUE 

EPS * SORT ( E . 0* C DUDX**S+0 . 5* C ( DUDY+DVDX )**S) +DVDy**S ) ) 

IF(EPS. LT . 0 . on EPS “ 0.01 
EPSN(KELEM,KIGS)« EPS 

VISC(KIGS) * ( 1 . 0+BSTAR*EPS**< 1 . 0/EN) )/t t • 73S*EPS) 

IFC lAXSY.EQ. 1 )DAREA»DAREA*RADUS 


C 

C««* 


PUT HEAT SEN, TERM i/or BODY FORCES INTO LOCAL RHS VECTOR. 



DO SO INOOP-1 , NNODP 

ERH5TCIN0DP.5 « ERHST < I NODP ) +SHAPP { I NODP 5 *EPS*PECLET*DAREA 
! *EPS*V] SCfKISS) 

60 CONTINUE 

DO ?0 ICON1=1 ,4 
DO 90 ICON£-1 , £ 

IROWU=( ICON! “I 3*7+4*ICONS-3 
IR0yV=lR0yU+3-IC0N2 
IROWT = IR0yU+4~IC0NE 
IN0DP=£*C1C0N1-1 )+ICON2 
SHAPI=SHAPP( INODP) 

CARXI=CARTP< 1 , INODP) 

CARf I=CARTPC£, INODP ) 

IFfICONE EQ.!) I ROWP* I ROyU+ 1 
DO 80 JC0N!*1 ,4 
JCOLP~f JCON1 -t )*?+£ 

KSALI = {KIGS-1 )*NNODL+JCONt 
C 

PUT PRESSURE TERM IN MOMENTUM EQUATIONS. 

C 

FLUMXC IROWU , JCOLP ) =FLUMX C I ROUU , JCOLP ) -CARXI *SHAl.G t KGAL I ) 

! *DAREA 

FLUMXi IROyV, JCOLP)*FLUf1X( IROyV, JCOLP)-CARYI*SHALS(KGALI )* 

« DAREA 

DO 80 JCONS»1 

JCOLU=( JCON1 “! )*7+4*JC0NE“3 
JC0LV=»JC0LU + 3~JC0N£ 

JCOLT » JCOLU + 4-JCONS 
JN0DP=£*{ JCON1 -1 }+JCONS 
SHAPJ=SHAPP{ JNODP) 

CARXJ=CARTP{ 1 , JNODP) 

CARY J = CARTP C£, JNOOP ) 

C 

C*** PUT DIFFUSION AND CONVECTION TERMS IN MOMENTUM EQUATIONS. 

C 

DIFFU*{ CARXI*E . 0*CARX J+CAR Y I *CARY J ) • VI SC ( K I SS ) *DAREA 
CONVC»( (UVELY*CARXJ+VVELY*CARYJ)*( SHAPI*DAREA) I/SISHA 
FLUMX( IROWU, JCOLU)==FLUHX{ IROWU, JCOLUJ+DIFFU+CONVC 
DIFFU « CCARXI»CARXJ+a. 0*CARYI*CARYJ)*VISCCKIGS)*DAREA 
FLUMXt IROWV, JCOLV)=FLUMX( IROUV, JCOLVI+DIFFU+CONVC 
FLUMXC IROWU, JCOLV) = CAR Y I *CARX J* VI SCI KI6S > *DAREA 
! +FLUMX( IROWU, JCOLV) 

FLUHXUROWV, JCOLU) «FLUHX{ IROWV , JCOLU)+ CARXI*CARYJ*VISC{ KIGS) 
? *DAREA 
C 

C* FORM ENERGY EQUATION 

C 

C*** PUT CONDUCTION AND CONVECTION TERMS IN THE ENERGY EQUATION 
C 

ENCOND ® CCARXI»CARXJ+CARYI*CARYJ)»DAREA 
ENCONV » <UVELY*CARXJ+VVELY*CARYJ)*SHAPI*DAREA»I‘PECLET 
FLUMX( IROUT, JCOLT) = FLUMXC IROUT, JCOLT )+ENCOND+ENCONV 
IFCICONE.EQ.EIGO TO 70 
C 

C*** form continuity EQUATION. 

C 

FLUHX< IROWP , JCOLU)»FLUMX( IROWP, JCOLU)+SHAPL« ICON) )*DAREA*CARXJ 
FLUMXI IROUP, JCOLV)aFLUMX{ IROWP, JCOLV )+SHAPL{ I CONI ) *DAREA*CARYJ 
FLUHX( IROWP , JCOLP )*0 . 0 
FLUMXC IROWP, JCOLT )*0. 0 
70 CONTINUE 

80 CONTINUE 

80 CONTINUE 

too CONTINUE 



c 

DO 1 1 0 JNODP’==1 ,NNODP 
KINP'-INODP 

KPOIf^^lABStLNODSC IELEM,KINP) ) 

I TOTU- NADFnCKPOIN) 

IT01V=1TCTU+N0DFM{KP0IN>-S 

irOTT ^ nOTU + NO0FM(KPOIN)-t 

EQPH54 I TOTU)=EQRHS{ ITOTU)+ERHSU{ INODP ) 

EQPH:.UTOTV>=EORHS{ I TOTV ) +ERHSV { INODP) 

EOPHSUTOTT) = EQRHSCITOTT) + ERHST(INODP) 

110 CONTINUE 

DO ^cO 1=1, H80UEL 
JK = I 

IFtKELEH.NE IBEL(JK)) GOTO tSO 
KBEL ~ lEEL(JK) 

KSURF =■ ISURFCJK) 

tFUKSURF.LE 3 ) . OR . ( KSURF . EQ . 7 ) ) GOTO ISO 
199 KSIDE ISIDECJK) 

TYPEx- , KSURF, KBEL, KS IDE 
type*, 'CALLING 80UC0N' 

CALL BDUCON (KBEL, KSURF, KSIDE, FLU«X,EQRHS,VARBE, HELEN, 

1 MPOIN,NNODP,MTOTV, IITER) 

ISO CONTINUE 

1987 TYPE 1988 

1988 FORHATOX, 'MATRIX OVER') 

RETURN 

END 

SUBROUTINE DR I VES ( COORD , LNODS , HELEN , MPO I N , NELEM , NGAUS , 

! NNOOL , NNODP, lELEH, AREAW) 

C*‘i'**»l<*'**’*** ***'**’•■ ■**♦#*«***♦>* 

c- 

C- EXTERNAL SUBROUTINES 

C- DJACOB 

C- 

C*»****«»*^f4-*#s + i»t%i|!i»'********: 

DIMENSION COORD{HPOIN,S) ,LN0DS<MELEH,8) 

D I HENS I ON DERIV(S,8) ,DJAC1{S,S) ,DJACK(E,E) ,5HAPE(8) 
DIMENSION CARTP<E,8),AREAWC9) 

COMHON/ARC1/POSCP(3) , WEICP«3) 

COMMON/ AREA3/CARPC(E, 72 ) ,SHAPG( 78) , CARLG ( S , 36 ) , SHALG ( 36 ) 
C COMHON/AREAA/AREAWI 9) 

C 

C*** REWIND TAPE BEFORE WRITING ON SHAPE FUNCTIONS 
C 

C REWIND 10 

C 

C*** SET UP POSITIONS AND WEIGHTS FOR 3 POINT GUASS RULE 
C 

POSGPi U«0 . T7A596669E 
P05GP( E)«0 . 0 
P0SCP(3)--P0SGP( 1 ) 

WCIGP( 1 )»0 .5555555556 
WE I GP ( E ) * 0 . 8888888889 
WEIGP? 3)«WeiCP( I ) 

C 

C*** CALCULATE SHAPE FUNCTIONS AND DERIVETIVES FOR ELEMENTS 
C 

LCAU5-0 

DO 50 IGAUS»1 ,NGAUS 
DO SO J6AUS»1 ,NGAUS 
LCAU8-LGAUS+1 
XE0IV«P0S6P(I6AU8) 

YEQlVxPOSCPI JGAUS) 

C 


c 

CALL .SHAPE8CDERiV,SHAPE,XEC!lV, YEQIV3 
C 

C*** SET UP JACOBIAN MATRIX AND INVERSE 
C 

CALL DJACOBCDERIV,DETJB,DJACI,DJACK, lELEM 
! ,MELEM,HPOIN,NNODP,LSAUS) 

C 

C*** CALCULATE GLOBAL DERIVETIVES AND AREA*GAUSS WEIGHTS 
C 

DO 10 IDIME=1,S 
DO 1 0 ■INODP = 1 ,NNODP 
CARTPC IDIME. INODP)=0 . 0 
DO 1 C JD3ME=1 ,£ 

CARTP( IDIME, INODP)=CARTP( IDIME, I NODP ) +D JAC 1 1 I D I ME , JD IME ) 
! *DERIVCJDiriE , INODP 3 
10 CONTINUE 

AREAUILGAUS )»DETJB*WEIGP( I GAUS ) *UE I GP C JGAUS 3 
C 

C»** PUT SHAPE FUNCTIONS AND DERIVETIVES IN ELEMENT MATRIX 
C 

DO 20 INODP=1 , NNODP 
KAPA= (LGAUS-1 3 *NNODP+ I NODP 
SHAPGI KARA )= SHAPE! INODP3 
DO £0 !D3ME*1 ,£ 

CARPG! I D IME, KARA )=CARTP( IDIME, INODP) 

£0 CONTINUE 

C 

C*** USE GAUSS POSITIONS TO CALCULATE LOCAL VALUES 
C 

CALL SHAPEA! DERI V, SHAPE, XEQIV,YEGIIV 3 
C 

C*»* CALCULATE GLOBAL DERIVATIVES FOR LINEAR FUNCTIONS 

C 

DO 30 IDIME»1 ,£ 

DO 30 INODL=1 ,NNODL 
CARTP< IDIME, INODLI«0. 0 
DO 30 JDIME=1 , £ 

CARTPC IDIME, INODL ) =CARTP ( IDIME, INODL3+ 

! DJACK IDIME, JDIHE3*DERIV( JDIHE, INODL3 
30 CONTINUE 

C 

C*** PUT SHAPE FUNCTIONS AND DERIVATIVES IN ELEMENT MATRIX 

C 

DO 40 INODL^I ,NNODL 
KGALI=<LGAUS-1 3 *NNODL+ I NODL 
SHALSCK6ALI }=SHAPE( INODL3 
DO 40 1D1HE=1 , £ 

CARLGC IDIME, KGAL I 3 “CARTPC IDIME, INODL3 
40 CONTINUE 

50 CONTINUE 

RETURN 
END 

SUBROUTINE SHAPES C DERI V, SHAPE , XEQIV , YEQIV 3 
DIMENSION DERIVCS,83 ,SHAPE(8) 

C 

C*«* PARABOLIC ELEMENT ANTICLOCKWISE NOD NUMBERING 

C 

X“XEQIV 

Y“YEQIV 

XY“X*Y 

XX»X*X 

YY*Y»»Y 

XXY*XX*Y 



X£=X-^£, 

Y£=Y*2 . 

XY2=XY*S 

SHAPE! n- (“t . +XY+XX+YY-XXY~XYY)* . 25 
SHAPE (£) = ( 1 ,-Y“XXHHXXY)*.5 
SHAPE ! 3 )-C -I . -XY+XX+YY-XXY+XYY J# . 25 
SHAPEC43 = {1 . +X-yY~XYY )* . 5 
SHAPE! 5 5 = ( “t . +XY>^XX+YY+XXY + XYY)* . 25 
SHAPE! 6 )= ( I ,+Y"XX-XXY)*.5 
SHAPE !7)= { -! . -XY+XX+YY+XXY-XYYJ* . 25 
SHAPE?85=(? .-X-YY+XYY)*.5 
DERIV! I . 1 )=•! Y+X2~XY2~YY)* . 25 
DERIVC t , £ )=-X+XY 
DEPI V( 1 , 3} = { ~Y+X2-XY2 + YY 3 * . 25 
OERIVCI ,4)={1 .-YY)*.S 
DERIV! t , 5?={ Y+X2+XYE+YY3* . 25 
DERIV! 1 ,6)*-X-XY 
DESIV! 3 ,?}=(-Y+X2+XY2-YY3*.25 
DERIV£1 ,S)=£-1 .+YY)*.5 
DERIV(2, 1 }*( X+V2-XX-XYE3* . 25 
C 

DERIV££,£3=£-1 .+XX)*.5 

DERIV(2,33»(-X+Y£-XX+XY23* .25 

OERIV(£,43=-Y-XY 

DERI VC 2, 5 3= (X + Y2+XX+XYE)* . 25 

DERIVtE, S 3 = C 1 .-XX3=»=.S 

DERIV(2,73=<-X+YE+XX“XY23*.E5 

DERIVtE, 8)=-Y+XY 

RETURN 

END 

|ri|e4i^Nle3^$2|(^#s^$:|;9#:9ica|i3fcs|cs|c^3|c3|i3|;#:|e9ies|eaiN:4cai!:|ca|ts|E:|!:|C)|( 

SUBROUTINE D JACOB ( DER I V , DET JB , D JAC I , DJACK , lELEH 
• ,MELEI1,f1POIN,NNODP,LGAUS) 

DIMENSION DERIV(E,83 , D JAC I ! 2 , 2 3 , D JACK < 2 , 2 3 

CO!1MON/ARCe/LNODS( 1 05,8 3 , COORD £ 366,2 > , NADFMt 366 3 , NODF«{ 366 3 
C 

C*** SET UP TEMPORARY MATRIX TO ALLOW THE JACOBIAN TO BE FORMED 
C 

DO 20 IDIHE=1 , 2 
DO £0 JDIME=1 ,2 
TEMPY=0 . £3 

DO 10 1N0DP»1 ,NNOOP 
KPOIN»IABS{LNODS( lELEM, INODP3 3 
KAG»CLGAUS-1 3 *NN0DP+1N0DP 

TEMPY=TEMPY+OERIV! IDIHE, 1N0DP3*C00RDCKP0IN, JDIME3 
10 CONTINUE 

DJACK( IDIME, JDIHE3*TEMPY 
20 CONTINUE 

DETJB«DJACKC 1 , 1 3 *0 JACK ( 2 , 2 3 -D JACK < 1 ,23*DJACK(2, 1 3 
C20S0 FORMAT(SX, 'DETERMENT VALUE" ' , FI 5 . 8 ) 

C 

C*** CHECK FOR NEGATIVE OR ZERO DETERMINENT 
C 

IF£DETJB330,30,40 
30 CONTINUE 

WRITE{S,20003 IELEH,L6AUS 
STOP 

40 CONTINUE 

C 

C*** INVERT TEMPORARY MATRIX TO FORM JACOBIAN 

C 

DJACIC 1 , ! 3»DJACKC2,23/DETJB 
DJACl C 2,23 "D JACK! 1 , 1 3/DETJB 
DJACI ( 1 ,2 3—DJACKC 1 ,2 3/DETJB 


y L- I . ‘ i * * W L 1 I B 

£000 F0RMArt//37H NON POSITIVE DETERMINENT FOR ELEMENT, EI 4 J 
RETURN 
END 

Kl- 4 **. it »*<!•* -s *******.♦*!<!*#*♦. 

SUBROUTINE SHAPE4 { DER I V , SHAPE , XEQ I V , YEQI V ) 

DIMENSION DERIVCa, 4) , SHAPEtA) 

C 

C«>it«««ti>* LINEAR ELEMENT ANTICLOCKWISE NODNUMBERINS 
C 

X=XEQIV 

Y=YEQIV 

XY=X*Y 

SHAPE { 1 )=< 1 . -X-Y+XY)* . £5 
SHAPE{£J*n . +X“Y-XY}* , 25 
SHAPEC3)-{ t .+X+y+XY)*.25 
SHAPED 4 ) = ( 1 .~X+Y-XYJ*.25 
DERIV( 1 , 1 }=(-! .+Y)*.25 
DERIVC 1 , a )*( 1 . -Y)* .25 
DERIV( 1 ,3)={ 1 .+YJ*.25 
DERIVI1 ,4J«C-1 .-Y)*.a5 
DERIV42, 1 ) = (-I .+X}*.Z5 
DERIV(2,2)=C-t .-Xi*.ZS 
DERIV{E,3)»( t .+XI*.25 
DERIV(2,4>=( 1 .-X)*.25 
RETURN 
END 

SUBROUTINE BOUCONIKBEL , KSURF , KSl DE , FLUHX , EQRHS , VARBE , 

1 MELEM,f1POIN,NNOBP,MTOTV, IITER) 

COnHON/ARC1/POSSP( 3> ,WEIGP(3) 

COMMON/ARCE/LNODSC 1 0S,8) , COORD (36 6, 2) ,NADFMC366 ) ,N0DFMC366) 
C0MM0N/ARC12/HSTAR,TFSTAR,FRIFAC,P1 ,EN, COSA,SINA, CONRAT 
COMMON/ARC1 3/USTAR , PECLET , SIGMA, RHO , MBOUEL 
COMMON/ARC30/ANGLNC366 I 

DIMENSION EGRHS(HTOTV) , VARBE ( HTOTV J , FLUHX ( 28 , 2S ) 

DIMENSION M(3 J ,N0DEC3> ,DERX(2,8) ,SHAP1 (8) ,SHAPN(8) 

DIMENSION RHSUCS) ,RHSV(5) ,RHST{83 , CART (2, 8) 

DIMENSION CJACKC2,£3 ,CJACI{2,2J ,C0SN{3I,SINN{3J 
C 

C** IDENTIFYING LOCAL NODE NOS. 

1F5KSURF.LE.3) SOTO 88 
H{ 1 3= \ 

MCE3 = 5 
M(3) = 8 

NODE( 1 3 = E*KSIDE - 1 
NODE(E3 = NODE! 1 3 + 1 
NODE(33 = NODE(£) + 1 
IF(NODE(33 .6T.85 NODEO) = NODE(33-8 
DO 145 I = 1,8 
RHSUC I 3 = 0.0 
RHSV(I3 = 0.0 
RHSTCn =0.0 
145 CONTINUE 

C — 

C GAUSS POINT EVALUATION STARTS 

C 

DO 150 IGAUS = 1,3 

GOTO (10,20,30,403 , KSIDE 



C CALCULATING LENGTH OF THE ELEMENT AT EACH GAUSS PT. 



^10 CONTINUE 

XEQIV = POSGP(IGAUS) 

YEQIV =-1.0 
I TEMP = 1 



aO CONTINUE 

XEQIi; *1.0 

YEQIV * POSGPdSAUSJ 

I TEMP * £ 

SOTO 50 

30 CONTINUE 

XEQIV = POSGPdGAUS) 

YEQIV =1.0 
ITEMP = 1 
SOTO SO 

40 CONTINUE 

XEQIV =-1.0 

YEQIV * POSSPiIGAUS) 

ITEMP = £ 

50 CONTINUE 

CALL SHAPE8CDERX, SHAP1 , XEQ I V , YEQ I V ) 

CALL DJACOB1 DERX. DETJB, CJACI , CJACK , KBEL , MELEH 
1 ,MPOIN,NWODP, IGAUSJ 

TEMPY = SQRT1CJACK( ITEMP, 1 >**£+CJACK( ITEMP, £}**£) 

SLETH = TEMPY««UEIGP( IGAUSJ 

C 

C CALCULATING VISCOSITY 4 STRAIN RATE AT GAUSS POINT FOR 5th SURF. ONLY. 

"" IFIKSURF . NE .S) GOTO 15 

DO 51 IDIME * 1 , S 
DO 51 INODP = I ,NNODP 
CART{ IDIME, INODP )«0 . 0 
DO 51 JDIME * 1,E 

CART I IDIME. INODP >«C ART ( IDIME, I NODP 1 +DERX ( JD I ME , I NODP ) 

1 *CJACI CIDIME, JDIME) 

51 CONTINUE 
UVELY =0.0 
VVELY =0,0 
DUDX =00 
DUDY =0.0 
DVDX =00 
DVDY =0.0 

DO 60 I = 1 , 8 

KPOIN = IABS(LNODS(KBEL, I ) ) 

ITOTU »NADF(1(KPOIN) 

ITOTV » ITOTU + NODFMIKPOIN) - £ 

SNAP * SHAP1 C I ) 

UVELY » UVELY + VARBS f I TOTU ) *SHAP 
VVELY « VVELY + VARSat ITOTV >*SHAP 
CARXI = CART! 1,1) 

CARYI = CART(a,I) 

UVEL = VARB£< ITOTU) 

VVEL = VARBaC ITOTV) 

DUDX = DUDX + CARX1*UVEL 
DUDY = DUDY + CARYI*UVEL 
DVDX = DVDX + CARXIoVVEL 
DVDY = DVDY + CARYI*VVEL 
60 CONTINUE 

EPS = SQRTCa. 0»(DUDX**a+0 .S*(DUDY+DVDX)**a+DV0Y»*a) ) 
IFIEPS.LT.O.OI ) EPS - 0.01 

Vise = (t .0+BSTAR#EPS*»<1 ,0/EN))/C1 .73a*EPS) 

VT « VVELYKCOSA + UVELY*SINA 



C EVALUATING GAUSS PT. INTEGRALS FOR LOCAL RHS VECTOR 

C 

DO 13 I = 1,3 
NOD » NODE(I) 

SHAP = SHAP! (NOD) 

FACTOR » FRIFAC*VI8C*EPS 

RHSU(NOD) « RHSU(NOD) “ SHAP*FACTOR*SI NA*SLETH 


RHSTfNOD'- * RHSTCNOD) + SHftP*PeCLET*FACTOR*VT*CONRAT*SLETH 
SHEAR * SHAF + PECLET*FAC.TOR*CONRAT*SLETH 
13 CONTINUE 

C SKIPPINS CONVECTIVE BC. FOR 5th SURFACE 

C 

IF(KSURF. EQ.S) GOTO 150 
15 CONTINUE 

C INCORPORATING APPROPRIATE CONVECTIVE HT . TR , BC. 

C 

DO 75 II = 1,3 

INODP = £1«{KSIDE-1 ) +I I 

IROWU = ?>(«{KSIDE-1 ) + Mdl) 

IFC INODP . GT . 8 J INODP = INODP - 8 
IF{ IROUU.GT.ES) IROWU a IROWU - 28 
KPOIN a lABSILNODSCKBEL, INODP) ) 

IROWT a IROWU + NODFHtKPOIN) -1 
SHAPI a shAPI ( INODP ) 

RHSTt INODP) a RHST( INODP) + SHAP 1 *HSTAR*TFSTAR*SLETH 

DO 95 JCONl a 1,4 

JCOLP a { JCONI-t >*7 + 2 

DO 9S JC0N2 a 1,2 

JNODP a ?JC0N1~1)*2 + JC0N2 

JCOLU a UCON1-n*7 + JC0N2*4-3 

JCOLT a UCON1-n*7 + JC0N£*3 + 1 

SNAP J a SHAPI ( JNODP ) 

TERN a HSTARtSHAPI*SHAPJ*SLETH 
FLUHXUROWT , JCOLT ) a FLUMX ( 1 ROWT , JCOLT ) +TERH 
95 CONTINUE 

TS CONTINUE 

ISO CONTINUE 

£6 CONTINUE 

C 

C INCORPORATING NORMAL VELOCITY BOUNDARY CONDITION 

C 

ANGLNC8£> = 90.0 
DO 80 II a 1,3 

C 

C EVALUATING DIRECTION OF NORMAL AT EACH NODE 

C — 

GOTO< 1 00, 1 1 0, 1 20, 130 } ,KSIDE 
too CONTINUE 

XEQIV a FLOATCin - 2.0 
YEQIV a -1 . 0 
ITEMP a 1 
RTEMP a +1 . 0 
GOTO 140 
110 CONTINUE 

XEQIV =1.0 

YEQIV a FLOATC II )~2. 0 

ITEMP a 2 

RTEMP =1.0 

GOTO 140 ' 

120 CONTINUE 

XEQIV a -FLOAT{II)+2.0 
YEQIV =1.0 
ITEMP a 1 
RTEMP a -1 .0 
GOTO 140 
130 CONTINUE 

XEQIV =“1.0 

YEQIV = -FLOAT( II )+£. 0 

ITEMP = 2 

RTEMP =-1.0 


call SHAF-E8CDERX, SHaPN,XEQIV, YEQIVI 

C EVALUATINGt^V^ff AT EACH NODE 

IHODP = £*?KSIDE-1 ) +11 
IROWU = 7*<KSIDE-1} + HUD 
IFdNODP.GT 8> INODP » INODP - 8 
IFCIROyy GT.£8) IROyU = IROyU - 28 
KPOIN = lABSCLNODStKBEL, INODP) ) 

IROWV - IROWU + NODFIKKPOIN) - £ 

IROyi = IROWV +1 

IFlKSURF.EQ.S) ANSLNCSfi) = £0.0 
ANGLE = ANGLN(KPOIN) 

ANGLE = ANGLE=)‘PI/180 . 0 
COSLX = RTEHP>isCOSCANGLE) 

COSNUI) = COSLX 

IFCABStCOSNC in ) .LT. 1 , OE-4) COSNdl) « 0.0 
COSLY * -RTEMP>^SIN{ ANGLE) 

SINNCin = COSLY 

IFt ABSCSINNt I n ) .LT . t . OE-4) SINNCII) = 0.0 
DO 90 JC0N1 = 1,4 
JCOLP = (JC0N1“U*7 + £ 

DO 90 JC0N2 =1,2 

JNODP = {JCON1-1)*£ + JCON£ 

JCOLU = CJCON1-1)*? + JC0N£*4-3 
JCOLV = JCOLU +3 - JCON£ 

JCOLT = JCOLU +4 - JCON£ 

SHAPPJ = SHAPNCJNODP) 

C 

C DETERNININS APPROPRIATE MOMENTUM EQN. FOR NORMAL VEL . BC . 


c 

1 

1 

1 

1 

1 

1 

1 

1 

i 

1 

1 

1 

f 

! 

i 

f 

1 

1 

I 

f 

1 

1 


I 

1 

I 

1 

1 

1 

1 

1 

1 

1 

1 

t 

1 

1 

1 

i 

1 

1 

i 

1 

1 

1 

! 

1 

1 

i 

I 

f 

I 


IFCABSCCOSNdl) ) . GE 

. ABSCSINNdl) ) ) GOTO 55 


FLUMX( IROWV, JCOLU) 

= COSNC 

I n*SHAPPJ 


FLUHXdROWV, JCOLV) 

= S1NN( 

1 1 )*SHAPPJ 


FLUMKdROWV, JCOLT) 

= 0.0 



IF( JCONE.EQ.S) GOTO 

56 



FLUHX( IROWV, JCOLP )“ 

0. 0 


58 

RHSVCINODP) =0.0 
GOTO 90 



5S 

CONTINUE 

FLUMXdROWU, JCOLU) 

= COSNC 

in*SHAPPJ 


FLUMX( IROWU, JCOLV) 

= SINNC 

II)*SHAPPJ 


FLUMXdROWU, JCOLT) 

= 0.0 



1F{ JCONE.EQ.E) GOTO 

67 



FLUMXdROWU. JCOLP) 

= 0.0 


67 

RHSUCINODP) = 0.0 



90 

CONTINUE 



8 0 

C 

CONTINUE 




C ASSEMBLING GAUSS PT . INTEGRALS IN GLOBAL RHS VECTOR 

C 

43 00 £00 II = 1,3 

NOD = NODEdl) 

KPOIN = 1ABS(LN0DS{KBEL,N0D) ) 

ITOTU = NADFMCKPOIN) 

ITOTV = ITOTU + NODFM(KPOIN) -2 
ITOTT « ITOTU + NODFMCKPOIN) -t 
EQRHSCITOTU) » EQRHS (ITOTU ) +RHSU( NOD ) 

EQRHSdTOTV) « EQRHSdTOTVI+RHSViNOD ) 

EORHSdTOTT) * EQRHSdTOTT>+RHST(NOD> 

£00 CONTINUE 

88 RETURN 

END 


